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y  Surface  tension  (dyne/cm,  1  dyne/cm=0. 00987  atm  pm) 
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INTRODUCTION 


In  vivo  gas  bubble  formation  and  growth  following  excessively  rapid  and  extensive 
decompression  is  considered  to  be  the  first  step  in  the  etiology  of  decompression 
sickness  (DCS)  in  man.1  Various  approaches  to  mitigating  the  incidence  and  severity  of 
DCS  consequently  entail  designing  exposures  in  which  individuals  must  be 
decompressed  to  limit  the  volumes  and  profusions  of  in  vivo  gas  bubbles.  Mathematical 
descriptions  of  gas  and  bubble  dynamics  in  tissue  are  central  to  such  approaches. 

Generally,  it  is  not  difficult  to  determine  the  temporal  course  of  bubble  growth  and 
resolution  in  a  well-stirred  solution,  especially  if  the  bubble  is  assumed  to  be  stationary.2 
In  contrast,  the  dynamics  of  gas  bubble  evolution  in  extravascular  tissue  -  evolution  that 
is  thought  to  give  rise  to  DCS  -  are  considerably  complicated  by  the  presence  of  blood 
flow.  A  detailed  mathematical  description  of  such  dynamics  involves  partial  differential 
equations,  the  solutions  of  which  are  far  too  complex  for  implementation  in  practical 
algorithms  for  mitigating  the  incidence  and  severity  of  DCS  in  humans.  Instead,  simpler 
models  have  been  developed  to  describe  the  underlying  processes  of  diffusion  and 
perfusion  in  terms  of  ordinary  differential  equations  (ODEs)  with  relatively  small 
numbers  of  parameters.  These  models  differ  in  the  way  bubble-tissue  gas  diffusion  and 
the  interactions  of  diffusion  processes  with  blood  flow  are  represented. 

This  report  reviews  a  series  of  ODE  models  of  gas  and  bubble  dynamics  in 
extravascular  tissue  that  were  examined  or  developed  in  earlier  work3'6  to  model  the 
evolution  of  in  vivo  gas  bubbles  thought  to  give  rise  to  DCS.  The  models  are  listed  in 
Tables  1  and  2.  All  except  the  two-region  (2R)  model  comprise  three  regions:  the  gas 
bubble,  a  diffusion  region  surrounding  the  bubble,  and  an  outer  tissue  region.  There  is 
no  outer  region  in  the  2R  model,  because  the  diffusion  region  is  theoretically  of  infinite 
extent.  In  the  three-region  models,  which  fall  into  two  classes  based  on  the 
characterization  of  gas  exchange  between  bubble  and  tissue,  the  outer  region  is 
considered  well-stirred,  implying  a  lack  of  gas  concentration  gradient  therein.  In  one 
class,  the  diffusion  region  around  the  bubble  is  a  thin  unperfused  layer  that  imposes  a 
barrier  to  gas  exchange  between  the  bubble  and  the  outer  well-stirred  tissue  region.  We 
refer  to  models  in  this  class  as  three-region  well-stirred  tissue  (3RWT)  models.  In  the 
other  class,  the  bubble  exchanges  gas  by  bulk  diffusion  only  with  tissue  in  the  middle 
diffusion  region,  which  is  relatively  large  and  includes  blood  flow:  The  outer  well-stirred 
region  does  not  participate  in  evolution  of  the  bubble.  We  show  that  the  2R  model  is  a 
special  case  of  this  model  class,  which  we  call  three-region  unstirred  tissue  (3RUT) 
models. 

We  present  the  models  in  the  logical  order  in  which  they  were  studied  or  developed, 
starting  from  a  three-region  well-stirred  tissue  model  with  a  constant  thickness  diffusion 
region  (3RWT-CT),  and  ending  with  a  three-region  unstirred  tissue  model  for  multiple 
bubbles  (3RUT-MB).  In  each  case,  the  presentation  includes  a  qualitative  description  of 
the  salient  features,  merits,  and  limitations  of  the  model.  The  Discussion  highlights  the 
merits  of  the  final  3RUT-MB  model. 
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Model  derivations  are  included  in  the  earlier  publications3'6  but  are  fractured  across  the 
different  articles  and  sketchy  because  of  space  limitations.  Consolidated  and  detailed 
derivations  of  the  model  equations  are  presented  in  Appendix  A,  with  inclusion  of  all 
mathematical  steps  involved.  The  presentation  is  aimed  at  providing  a  better 
understanding  of  the  differences  between  the  various  models,  differences  that  are 
obscured  by  the  structural  similarities  of  the  final  model  equations.  No  numerical  results 
are  presented,  as  detailed  simulation  results  for  each  of  the  models  may  be  found  in  the 
earlier  publications.  Appendix  B  presents  a  discussion  of  the  correspondence  between 
the  diffusion  and  perfusion  equations  in  which  it  is  shown  that  the  latter  is  simply  the 
integral  version  of  the  former,  which  underscores  the  need  to  include  perfusion  terms  in 
the  diffusion  equation  and  bubble-tissue  diffusion  terms  in  the  perfusion  equation.  A 
brief  overview  of  tissue  elastic  effects  on  bubble  pressure  is  presented  in  Appendix  C. 


MATHEMATICAL  MODELS  OF  GAS  BUBBLE  EVOLUTION 
Well-stirred  Tissue  Models 


I.  Three-Region  Well-Stirred  Tissue  Model  with  Constant  Thickness  Diffusion 
Region  (3RWT-CT) 

The  three-region  well-stirred  tissue  (3RWT)  model  consists  of  the  gas  bubble,  an 
unperfused  thin  diffusion  layer  surrounding  the  bubble,  and  an  outer  well-stirred  tissue 
region.  The  diffusion  layer  is  a  lumped  representation  that  accounts  for  all  gas  diffusion 
distributed  in  the  tissue.  The  gas  flux  into  and  out  of  the  bubble  depends  on  the 
diffusivity  of  this  layer.  The  gas  fluxes  at  both  boundaries  of  the  diffusion  layer  are  the 
same  at  any  time  during  bubble  growth,  because  of  uniform  diffusivity  and  lack  of  blood 
flow  in  the  layer.  Gas  concentration  gradients  exist  only  in  the  diffusion  layer,  because 
the  region  outside  the  layer  is  well-stirred. 

The  3RWT  model  with  a  constant  diffusion  layer  thickness  (3RWT-CT)  that  we 
described  in  the  first  paper  of  our  series3  is  a  modified  version  of  the  model  described 
by  Gernhardt.7  The  two  models  differ  in  two  important  aspects.  The  present  version 
accounts  for  gas  losses  or  gains  by  the  bubbles  in  computing  the  tissue  gas  tension. 
Secondly,  without  any  additional  parameter,  the  3RWT-CT  model  eliminates  the 
approximation  in  the  Gernhardt  model  that  the  diffusion  layer  thickness  h  is  much 
smaller  than  the  bubble  radius  n. 

We  demonstrated  the  significant  effect  of  bubble  growth  on  tissue  gas  tension  at  low 
tissue  volumes  in  a  previous  publication  (Fig.  2,  Ref.  3).  For  a  given  tissue  volume, 
there  is  a  limit  to  bubble  number  density  above  which  the  influence  of  the  bubbles  on 
tissue  gas  tension  cannot  be  ignored.  We  also  discussed  the  result  of  assuming  h  «  n 
in  the  same  article  (Fig.  4,  Ref.  3).  With  this  approximation,  h  cannot  be  3  pm,  as 
Gernhardt  assumed  in  applications  of  his  model,  if  bubble  growth  is  started  from  an 
initial  bubble  nucleus  radius  of  3  pm  or  less.  The  approximation  may  not  significantly 
affect  the  calculation  of  bubble  radius,  but  that  depends  on  the  particular 
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decompression  profile  used.  In  any  case,  the  approximation  is  not  necessary  and  is 
eliminated  in  our3RWT-CT  model. 

The  drawbacks  of  the  3RWT-CT  model  are  (i)  the  diffusion  layer  is  a  factitious  layer  with 
an  arbitrary  thickness;  (ii)  as  a  result  of  equal  gas  fluxes  at  the  inner  and  outer 
boundaries  of  the  diffusion  layer,  the  gas  content  of  the  layer  must  remain  unchanged 
as  the  ambient  pressure  changes  during  decompression,  which  causes  the  dissolved 
gas  in  the  layer  to  violate  Henry’s  law;  and  (iii)  zero  concentration  gradient  in  the  well- 
stirred  region  leads  to  an  abrupt  change  in  gas  flux  from  zero  just  outside  the  diffusion 
layer  to  a  finite  value  just  inside  the  layer. 

Despite  the  inconsistencies,  the  3RWT-CT  model  is  a  multiple  bubble  model  that  takes 
into  account  the  competition  between  diffusion  and  perfusion  processes  for  available 
gas.  As  demonstrated  in  our  original  consideration  of  this  model,3  the  increased  rates  of 
gas  depletion  that  occur  with  growth  of  bubbles  at  increased  bubble  number  densities 
can  severely  limit  the  size  attained  by  the  bubbles.  At  their  extreme  such  effects  cause 
the  bubble  gas  pressures  to  be  clamped  to  the  ambient  hydrostatic  pressure,  so  that 
during  subsequent  resolution  of  the  bubbles  gas  elimination  from  the  tissue  at  constant 
ambient  pressure  becomes  a  linear  function  of  time.  These  effects,  however,  are  not 
manifestations  of  direct  bubble-bubble  interactions,  which  are  not  accommodated  until 
we  consider  the  3RUT-MB  model  in  Section  V  below. 

An  abrupt  transition  of  gas  flux  at  an  arbitrary  distance  from  the  bubble  surface,  where 
there  is  not  necessarily  any  coinciding  physical  surface,  is  an  inherent  conceptual 
problem  with  all  well-stirred  tissue  models.  The  violation  of  Henry’s  law  in  the  diffusion 
layer  is  perhaps  a  more  disconcerting  problem,  but  one  that  can  be  rectified  by  allowing 
the  thickness  of  the  diffusion  layer  to  vary  during  bubble  growth. 

II.  Three-Region  Well-Stirred  Tissue  Model  with  Varying  Thickness  Diffusion 
Region  (3RWT-VTDD) 

In  this  version  of  the  3RWT  model,4  the  thickness  of  the  diffusion  region  surrounding  the 
gas  bubble  is  allowed  to  vary  by  postulating  a  difference  in  gas  diffusivity  between  an 
infinitesimally  thin  layer  at  the  bubble  surface  and  the  remainder  of  the  diffusion  region. 
We  refer  to  this  model  as  the  varying  thickness  differential  diffusivity  (VTDD)  version  of 
the  3RWT  model.  The  changes  in  diffusion  region  thickness  allow  dissolved  gas  content 
in  the  region  to  change  in  accord  with  Henry’s  law,  thus  eliminating  a  theoretical 
inconsistency  associated  with  the  3RWT-CT  model. 

Comparison  with  the  3RWT-CT  model  shows  the  behavior  of  the  3RWT-VTDD  model  to 
be  more  regular  and  consistent  across  physiological  ranges  of  the  model  parameters. 
Also,  the  3RWT-VTDD  model  predicts  longer  bubble  lifetimes  in  tissue  compartments 
with  shorter  half-times  than  those  obtained  using  the  3RWT-CT  model  (Figs.  5  and  6, 
Ref.  4). 


3 


The  differential  diffusivity  does  not  introduce  any  additional  parameter  in  the  3RWT- 
VTDD  model,  which  is  structurally  identical  to  the  3RWT-CT  model.  The  constant 
thickness  is  simply  replaced  by  a  variable  thickness  determined  by  the  ratio  of  diffusion 
region  diffusivity  to  bubble  surface  diffusivity.  But  additional  computation  is  required  to 
solve  a  cubic  polynomial  for  the  diffusion  region  thickness  at  every  step  of  bubble 
growth  and  resolution. 

In  short,  the  3RWT-VTDD  model  replaces  the  factitious  diffusion  layer  of  the  3RWT-CT 
model  by  an  equally  factitious  bubble  surface  diffusivity  that  helps  to  preserve  mass 
balance  in  the  diffusion  region  without  violating  Henry’s  law.  In  the  absence  of  any 
biophysical  mechanism  that  may  be  invoked  to  justify  its  existence,  differential  diffusivity 
of  gas  at  the  bubble  surface  is  only  a  convenient  mathematical  construct.  The  3RWT- 
VTDD  model  also  features  the  abrupt  transition  of  gas  flux  at  the  outer  boundary  of  the 
diffusion  region  that,  as  mentioned  earlier,  is  characteristic  of  all  well-stirred  tissue 
models.  A  more  serious  inconsistency  is  that  the  diffusion  layer  is  no  longer  a  thin  layer 
around  the  bubble,  but  can  expand  into  a  large  unperfused  diffusion  region  during 
bubble  growth,  depending  on  the  differential  diffusivity  ratio.  The  existence  of  such 
unperfused  diffusion  regions  around  bubbles  in  actual  tissue  becomes  ever  more 
untenable  as  the  sizes  of  these  regions  increase.  It  is  equally  disconcerting  that  when 
the  model  is  elaborated  to  accommodate  multiple  diffusible  gases,  the  diffusion  region 
thickness  and  volume  around  each  bubble  must  be  different  for  each  gas  due  to  the 
different  diffusivity  ratios,  tissue  gas  tensions,  and  bubble  pressures  of  each  gas. 


Unstirred  Tissue  Models 


The  inconsistencies  and  additional  computational  requirements  of  3RWT  models  forced 
us  to  turn  our  attention  to  unstirred  tissue  models  of  gas  bubble  dynamics  in  which  gas 
exchange  is  limited  by  bulk  diffusion  through  the  tissue.  These  models  explicitly  account 
for  the  effects  of  perfusion  on  gas  diffusion.  The  capillaries  are  considered  as  sinks  for 
gas  diffusing  out  of  the  bubble  or  as  sources  for  gas  diffusing  into  the  bubble.  Realistic 
representation  of  the  sinks  as  discrete  regions  distributed  throughout  the  tissue  entails 
excessive  computational  complexity.  Their  representation  is  simplified  by  assuming  that 
the  blood  flow  is  uniform  and  present  everywhere  in  the  diffusion  region.  This 
assumption  poses  no  inconsistency  as  we  can  calculate  the  requisite  blood  flow  per  unit 
volume  of  the  diffusion  region  based  on  the  total  blood  flow  to  the  region.  Also,  there  is 
no  requirement  for  the  diffusivity  of  the  gas  in  bulk  tissue  to  differ  from  the  diffusivity  of 
the  gas  at  the  bubble  surface.  Finally,  the  assumption  of  perfusion  uniformity  can  be 
relaxed  to  include  the  effects  of  perfusion  heterogeneity,  as  discussed  in  connection 
with  the  3RUT  models  in  sections  IV  and  V  below. 

The  gas  lost  in  the  sinks  is  calculated  from  the  difference  between  a  sink  pressure  and 
the  prevailing  tissue  gas  tension.  The  sink  pressure  would  be  the  arterial  tension  of  the 
gas,  if  the  sinks  are  represented  as  distinctly  distributed  regions.  But  with  presumption 
of  a  continuous  distribution  of  blood  flow  in  the  diffusion  region,  we  define  the  sink 
pressure  as  the  spatial  average  gas  tension  Pd  in  the  region. 
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III.  Two-Region  Model  (2R) 

The  2R  model  is  an  unstirred  tissue  model  of  gas  bubble  dynamics  that  consists  of  only 
two  regions;  bubble  and  tissue.  It  was  originally  formulated  by  Van  Liew  and  Hlastala 8 
to  simulate  the  dissolution  of  subcutaneous  gas  pockets  in  rats  and  has  since  been 
used  to  simulate  the  dynamics  of  in  vivo  gas  bubbles  in  animals  and  humans.9'11  The 
entire  tissue  volume  participates  in  gas  diffusion  in  the  2R  model,  with  the  spatial 
average  gas  tension  Pd  in  this  volume  serving  as  the  sink  pressure.  It  is  shown  in 
Appendix  A  that  the  uniform  tissue  gas  tension  Pt  far  away  from  the  bubble,  which  is  the 
sink  pressure  assumed  by  Van  Liew  and  Hlastala  in  their  2R  model  formulation,  is  the 
same  as  this  average  gas  tension  Pd . 

In  an  earlier  article,3  we  pointed  out  the  limitations  of  the  2R  model.  The  assumption  of 
uniform  sink  pressure  along  with  the  boundary  condition  that  the  gas  concentration 
gradient,  and  hence  the  gas  flux,  vanish  at  the  periphery  requires  the  tissue  volume 
around  the  bubble  to  be  theoretically  infinite.  Therefore,  the  model  is  not  applicable  for 
evaluating  the  evolution  of  multiple  bubbles  or  for  evaluating  the  evolution  of  even  a 
single  bubble  in  a  finite  volume  of  tissue.  Also,  uniform  blood  flow  implies  a  uniform 
distribution  of  capillaries  in  the  unstirred  tissue,  which  ignores  the  possibility  that 
perfusion  can  be  heterogeneous. 

IV.  Three-Region  Unstirred  Tissue  Model  (3RUT) 

The  shortcomings  of  the  3RWT  model  stem  chiefly  from  a  non-zero  gas  concentration 
gradient  at  the  boundary  between  the  diffusion  region  and  the  outer  well-stirred  region 
around  a  bubble.  If  this  gradient  can  be  made  zero,  all  gas  exchange  is  strictly  between 
the  bubble  and  the  unstirred  diffusion  region  around  it,  with  no  role  for  gas  in  the  outer 
well-stirred  region.  This  is  accomplished  in  the  three-region  unstirred  tissue  (3RUT) 
model  by  introducing  deviations  to  the  uniform  sink  pressure,  deviations  that  in  turn 
cause  gas  tension  deviations  in  the  diffusion  region.5  The  model  exploits  the  fact  that 
only  the  gas  concentration  gradient  at  the  bubble  surface  is  needed  to  determine  the 
temporal  course  of  bubble  growth  and  resolution.  The  sink  pressure  deviations  are 
represented  such  that  they  do  not  affect  the  gas  flux  at  the  bubble  surface  and  also 
ensure  a  zero  gas  flux  at  the  outer  boundary  of  the  diffusion  region,  which  is  assumed 
to  be  of  constant  volume  during  the  lifetime  of  the  bubble. 

The  gas  flux  at  the  bubble  surface  in  the  3RUT  model  is  determined  by  the  spatial 
average  gas  tension  Pd  in  the  diffusion  region,  which  need  not  be  the  same  as  the  gas 
tension  Pt  in  the  surrounding  well-stirred  region.  As  the  diffusion  region  volume 
increases,  the  sink  pressure  deviations  decrease,  and  the  model  reduces  to  the  2R 
model  with  Pd  ->  Pt.  The  computations  require  specification  of  only  the  volume  of  the 
diffusion  region,  not  its  thickness.  This  volume  must  be  above  a  certain  minimum  value 
in  order  to  sustain  bubble  growth,  because  all  gas  exchange  takes  place  exclusively 
from  the  diffusion  region.  The  adequacy  of  the  diffusion  region  volume  must  be  checked 
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during  bubble  growth  and  resolution,  but  the  needed  additional  computations  are 
minimal. 

The  strength  of  the  3RUT  model  lies  in  the  representation  of  sink  pressure  that  allows 
the  diffusion  region  volume  to  be  finite  without  compromising  overall  computational 
simplicity.  More  importantly,  the  deviations  in  sink  pressure  can  account  for  complex 
heterogeneous  perfusion  effects  on  gas  bubble  dynamics.  However,  because  the 
diffusion  region  around  the  bubble  is  assumed  to  be  radially  symmetric,  the  outer 
boundary  of  the  diffusion  region  must  everywhere  adjoin  the  outermost  well-stirred 
tissue  region,  which  does  not  participate  in  bubble  evolution,  to  meet  the  requirement 
for  zero  flux  at  the  boundary.  As  a  result  of  this  constraint  and  the  assumption  that  the 
volume  of  the  diffusion  region  around  each  bubble  is  constant,  each  bubble  grows  and 
resolves  in  its  own  independent  diffusion  region,  which  accounts  for  the  finiteness  of 
tissue  size  but  precludes  consideration  of  the  increased  rates  of  gas  depletion  that 
occur  with  growth  of  bubbles  at  increased  bubble  number  densities  or  direct  bubble- 
bubble  interactions.  Like  the  2R  model,  the  3RUT  model  is  applicable  to  the  evolution  of 
only  a  single  bubble  in  a  tissue. 

V.  Three-Region  Unstirred  Tissue,  Multiple  Bubble  Model  (3RUT-MB) 

The  3RUT  multiple  bubble  (3RUT-MB)  model 6  further  exploits  the  fact  that  only  the  gas 
concentration  gradient  at  the  bubble  surface  is  required  to  determine  bubble  growth  and 
resolution.  The  concept  of  sink  pressure  deviations  is  generalized,  assuming  that  only 
the  components  of  the  gas  concentration  gradient  normal  to  the  bubble  surface 
contribute  to  the  average  gas  flux  into  and  out  of  the  bubble.  This  implies  that  (i)  gas 
diffusion  need  not  be  spherically  symmetric  around  the  gas  bubble;  (ii)  the  diffusion 
region  around  the  bubble  may  be  of  any  shape;  and  (iii)  the  diffusion  region  boundaries, 
and  hence  its  volume,  may  vary  during  the  course  of  bubble  evolution.  Because  the 
diffusion  region  around  a  given  bubble  can  be  spherically  asymmetric,  the  requirement 
for  zero  flux  at  the  outer  boundary  can  be  met  without  requiring  the  outer  boundary  to 
everywhere  adjoin  the  outermost  well-stirred  tissue  region,  as  in  the  3RUT  model.  The 
diffusion  region  of  one  bubble  can  consequently  abut  the  diffusion  region  of  one  or  more 
adjoining  bubbles,  allowing  direct  bubble-bubble  interactions  that  cannot  be 
accommodated  in  the  3RUT  model.  The  diffusion  region  gas  tension  is  determined 
based  on  an  average  blood  flow  assuming  that  the  net  gas  content  carried  by  blood  into 
and  out  of  the  tissue  by  the  heterogeneous  perfusion  components  is  zero. 

While  the  diffusion  region  volume  around  any  given  bubble  can  vary  during  the  lifetime 
of  the  bubble  in  the  3RUT-MB  model,  the  total  tissue  volume  exclusive  of  bubble 
volumes  does  not  vary.  At  any  time  during  bubble  growth  or  resolution,  this  total  volume 
is  equal  to  the  sum  of  the  diffusion  region  volumes  of  all  bubbles  plus  the  volumes  of 
any  bubble-free  regions  within  the  tissue.  This  total  volume  places  an  upper  limit  on  the 
bubble  number  density  that  a  given  volume  of  tissue  can  accommodate,  analogous  to 
the  lower  limit  on  the  diffusion  region  volume  in  the  3RUT  model.  This  density  limit  must 
be  checked  during  bubble  growth  and  resolution,  but  again  the  needed  additional 
computations  are  minimal. 
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The  bubble-bubble  interactions  accommodated  by  the  3RUT-MB  model  are  different 
from  and  act  in  addition  to  the  increased  rates  of  gas  depletion  that  occur  with  growth  of 
bubbles  at  increased  bubble  number  densities.3  The  latter  effects  are  accommodated  by 
this  model  because  as  the  diffusion  region  volumes  of  the  bubbles  can  change 
continuously  during  bubble  evolution,  the  bubbles  share  the  same  average  diffusion 
region  gas  tension  which  governs  the  rate  of  bubble  growth  or  resolution.  The  additional 
bubble-bubble  interactions  accommodated  by  the  model  include  proximity-driven 
reductions  of  the  diffusion  region  volumes  of  adjoining  bubbles.  Such  reductions  can 
further  impair  the  growth  of  bubbles  whose  growth  at  high  bubble  number  densities  is 
already  impaired  by  the  accelerated  depletion  of  available  gas.  The  reductions  of 
diffusion  region  volumes  that  cause  these  impairments  are  the  3RUT-MB  model 
representations  of  the  bubble-bubble  interactions  modeled  in  detail  by  Jiang,  et  al.,12  in 
a  two-bubble  system  without  perfusion  or  mass  balance. 


DISCUSSION 

A  principal  objective  of  the  model  examination  and  development  effort  summarized  here 
was  to  identify  a  system  of  equations  for  modeling  the  evolution  of  gas  bubbles  in 
tissue,  equations  that  remain  sufficiently  simple  to  be  applicable  in  algorithms  for 
mitigating  risks  of  DCS  in  diving,  altitude,  or  flying-after-diving  profiles,  and  that  are 
based  on  a  realistic  conceptualization  of  the  bubble-tissue  system.  We  have  always 
presumed  that  an  algorithm  based  on  sound  biophysical  principles  will  not  only  allow 
correlation  of  DCS  incidence  and  time  of  occurrence  data  from  the  widest  possible 
variety  of  decompression  profiles,  but  also  allow  the  risks  of  DCS  in  profiles  different 
from  those  included  in  the  correlation  to  be  evaluated  with  higher  confidence. 

Tables  1  and  2  summarize  the  features  and  limitations  of  the  various  ODE  models  of 
diffusion-limited  gas  bubble  dynamics  considered  in  this  work  and  show  the 
corresponding  model  equations  for  computing  the  evolution  of  gas  bubble  volume  and 
tissue  gas  tension.  The  equations  for  the  different  models  are  strikingly  similar,  despite 
the  very  different  conceptualizations  of  the  bubble-tissue  system  in  the  different  models. 
This  is  partly  because  all  equations  are  based  on  a  spherically  symmetric  solution  of  the 
diffusion  equation  under  the  quasi-static  approximation,  and  partly  by  design  following 
from  appropriate  formulation  of  the  simplifying  assumptions. 

In  all  models,  the  expression  for  the  gas  tension  gradient  at  the  bubble  surface  includes 
a  1/n  term,  which  embodies  the  divergence  of  the  gradient,  modified  by  addition  of 
another  term  that  embodies  all  model-specific  features  of  gas  diffusion  between  bubble 
and  tissue.  This  added  term  is  a  fixed  parameter  in  the  2R,  3RWT-CT,  and  3RUT 
models  if  the  tissue  perfusion  rate  remains  constant,  but  generally  varies  with  time  in 
the  remaining  3RWT-VTDD  and  3RUT-MB  models.  Such  time  dependence  is  explicitly 
defined  in  the  3RWT-VTDD  model. 
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All  models  except  the  single-bubble  2R  and  3RUT  models  are  applicable  to  modeling 
the  evolution  of  in  vivo  extravascular  bubbles.  Of  the  applicable  models,  the  3RWT-CT 
and  3RUT-MB  models  are  structurally  identical  if  the  term  added  to  1/n  in  the  dr/dt 
equation  for  each  of  the  two  models  is  assumed  constant.  Under  this  assumption,  either 
of  these  models  will  prove  equally  able  to  correlate  data  for  the  incidence  and  time  of 
occurrence  of  DCS  in  diving,  altitude,  or  flying-after-diving  profiles  when  used  in  the 
same  fashion  to  track  the  accumulation  of  DCS  risk  in  such  profiles.  However,  the 
interpretation  and  acceptable  ranges  of  the  model  parameter  values  remain 
unconfounded  by  violation  of  physical  laws  only  in  the  3RUT-MB  model. 

The  assumptions  underlying  the  3RUT-MB  model  are  not  only  more  theoretically 
consistent  than  those  underlying  the  3RWT-CT  model,  but  are  also  more  conceptually 
realistic.  The  3RUT-MB  model  is  based  on  the  assumption  that  extravascular  bubbles 
exchange  gas  by  diffusion  through  an  unstirred  and  heterogeneously  perfused  tissue,  in 
accord  with  the  intrinsically  heterogeneous  nature  of  tissue  perfusion.  The  parameter 
added  to  1/n  in  the  expression  for  the  gas  tension  gradient  at  the  bubble  surface 
embodies  not  only  the  influences  of  the  bulk  diffusivity  of  gas  in  the  tissue  and  the 
blood-tissue  gas  exchange  time  constant,  but  also  the  effects  of  perfusion  heterogeneity 
and  interactions  between  multiple  bubbles.  The  widest  possible  scope  of  such  effects  is 
accommodated  as  each  bubble  is  not  considered  to  be  at  the  center  of  its  own  diffusion 
space  in  the  tissue  with  gas  diffusion  in  only  the  radial  direction  from  the  bubble  center, 
but  is  instead  considered  to  be  in  an  environment  in  which  gas  can  diffuse  in  arbitrary 
directions  around  the  bubble.  Thus,  while  this  parameter  is  constrained  in  the  3RWT-CT 
model  by  its  identification  with  the  thickness  of  the  diffusion  layer  around  the  bubble,  the 
corresponding  A  parameter  in  the  3RUT-MB  model  can  be  specified  without  such 
constraint  and  can  validly  assume  a  much  wider  range  of  values. 

If  constancy  of  the  term  added  to  1/n  in  the  dr/dt  equation  for  the  3RWT-CT  and  3RUT- 
MB  models  is  relaxed,  variations  of  this  parameter  in  3RWT-CT  model  remain 
constrained  by  its  identification  with  the  thickness  of  the  diffusion  layer  around  the 
bubble:  The  parameter  can  only  vary  within  reasonable  ranges  of  diffusion  layer 
thicknesses.  While  this  parameter  is  always  time-dependent  in  the  remaining  applicable 
3RWT-VTDD  model,  its  variations  are  similarly  constrained  as  it  is  an  explicit  function  of 
the  size  and  content  of  the  bubble  as  well  as  the  tissue  gas  tension.  Again,  the 
corresponding  A  parameter  in  the  3RUT-MB  model  can  vary  without  such  constraints. 
Time-dependent  variations  in  this  parameter  afford  the  added  flexibility  to  accommodate 
the  influences  of  variations  in  perfusion  heterogeneity  and  bubble-bubble  interactions 
that  may  be  expected  as  the  bubble  sizes,  the  work  and  thermal  status  of  an  exposed 
individual,  and  the  number  of  bubbles  in  the  tissue  change  during  a 
compression/decompression  profile.  Accommodation  of  such  variations  is  beyond  the 
scope  of  the  other  models. 

The  value  of  A  in  the  3RUT-MB  model  cannot  be  specified  a  priori  without  a  detailed 
description  of  the  diffusion  field  around  each  bubble,  but  is  readily  determined  by  fitting 
the  model  to  empirical  data  assuming  that  A  is  the  same  for  all  bubbles.  Similarly,  time 
dependence  of  A  can  be  accommodated  by  incorporation  of  empirical  relationships 
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between  A  and  the  compartmental  blood-tissue  gas  exchange  time  constant,  bubble 
size,  and  bubble  number  density.  With  such  accommodation,  the  environment 
governing  bubble  evolution  can  change  within  a  compartment  of  given  blood-tissue  gas 
exchange  halftime  with  properties  of  the  pressure/time/breathing  gas  profile,  greatly 
increasing  overall  model  flexibility. 


CONCLUSION 

The  3RUT-MB  model  is  preferred  over  the  alternative  models  for  its  enhanced 
conceptual  realism  and  flexibility.  Time  dependence  of  the  A  parameter  should  be 
empirically  accommodated  to  account  for  variations  in  perfusion  heterogeneity  and 
bubble-bubble  interactions  that  may  be  expected  as  the  bubble  sizes,  the  work  and 
thermal  status  of  an  exposed  individual,  and  the  number  of  bubbles  in  the  tissue  change 
during  a  compression/decompression  profile. 
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Table  1 .  Models  of  bubble  evolution  in  well-stirred  extravascular  tissue 


3RWT-CT  -  Constant  Diffusion  Region  Thickness 


Description 

Gas  bubble  is  immersed  in  a  well-stirred 
tissue  compartment  of  finite  volume  and  is 
immediately  surrounded  by  an  unperfused 
boundary  layer  of  constant  thickness 
through  which  diffusion-limited  gas 
exchange  between  bubble  and  well-stirred 
tissue  occurs. 


Equations 

Rate  of  change  of  radius  of  each  identical¬ 
sized  bubble: 


dr, 

dt 


a'tDs 


1  1 

— +  — 

h  r, 


(pt-pi)-7d7f 

3  dt 


P  —  P  +  ^  ^  Mr  ^ 

*amb  Mdg 

3r,  3 


Features 

•  Spherically  symmetric  diffusion. 

•  Applicable  to  multiple  bubbles. 

•  Dissolved  gas  in  diffusion  layer  violates 
Henry’s  law. 

•  Discontinuity  of  gas  flux  at  outer 
boundary  of  diffusion  layer. 

•  Small  diffusion  layer  volume  compared 
to  tissue  volume;  may  be  ignored  at  low 
bubble  number  densities. 

•  Necessary  to  account  for  diffusion  layer 
volume  in  determining  bubble-tissue  gas 
exchange  volume  at  higher  bubble 
number  densities. 


where  h  is  a  constant  diffusion  region 
thickness. 

Tissue  gas  tension  with  multiple  bubbles: 


dgt  _  Pa  ~Pt 

dt  x 

neglecting  diffusion 


atvt  sdt 


layer  volume. 
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Table  1 .  ( continued ) 


3RWT-VTDD  -  Varying  Diffusion  Region  Thickness 

Description 

Gas  bubble  is  immersed  in  a  well-stirred 
tissue  compartment  of  finite  volume  and  is 
immediately  surrounded  by  an  unperfused 
boundary  layer  of  varying  thickness 
through  which  diffusion-limited  gas 
exchange  between  bubble  and  well-stirred 
tissue  occurs. 

Features 

•  Spherically  symmetric  diffusion. 

•  Applicable  to  multiple  bubbles. 

•  Dissolved  gas  in  diffusion  layer  obeys 
Henry’s  law. 

•  Diffusion  layer  thickness  determined  by 
bubble  size  and  content,  the  ratio  p  of 
the  bulk  tissue  gas  diffusivity  to  the 
bubble  surface  gas  diffusivity,  and  the 
gas  tension  Pt  in  the  well-stirred  region. 

•  Discontinuity  of  gas  flux  at  both 
boundaries  of  diffusion  layer. 

•  Diffusion  layer  volume  comparable  to 
tissue  volume;  correction  may  be 
needed  even  at  low  bubble  number 
densities. 

Equations 

Rate  of  change  of  radius  of  each  identical¬ 
sized  bubble: 

dr,  |_h  l-J  3  dt 

dt  ~  P  -P  +4y  +  8ltMr3 

where  h  is  the  diffusion  region  thickness 
determined  by  solving  a  cubic  polynomial 
in  h. 

Tissue  gas  tension  with  multiple  bubbles: 

dPt_Pa-Pt  P  V  d  /p  y  j 

dt  t  atVt(l-v)  j—f  dt 

where  v  =  ^(Vd)k  /vt  is  the  ratio  of  total 

k=l  / 

diffusion  layer  volume  from  all  bubbles  to 
tissue  volume;  with  valid  bubble  number 
densities,  v  <  1 . 

The  diffusion  region  thickness  hk  for  each 
of  k=1,  ...,  N  bubbles  and  v  are  updated  at 
the  beginning  of  each  integration  step. 
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Table  2.  Models  of  bubble  evolution  in  unstirred  extravascular  tissue 


2R  -  Infinite  Tissue  Volume 


Description 

Gas  bubble  is  immersed  in  an  unstirred 
uniformly  perfused  tissue  compartment  of 
infinite  extent,  through  which  diffusion- 
limited  gas  exchange  between  bubble  and 
tissue  occurs. 

Features 

•  Spherically  symmetric  diffusion. 

•  Tissue  volume  is  theoretically  infinite, 
and  is  not  a  factor  in  determining  bubble 
growth.  Corollary: 

•  Applicable  only  to  a  single  bubble 
immersed  in  a  physiologically 
unrealistic  infinite  volume  of  tissue. 


Equations 

Rate  of  change  of  bubble  radius: 


dt 


a'tDs 


X  +  - 


(Pt-PiK 


r  *fPamb 


3  dt 


P  —  P  4y  Mr3 

*amb  Mdg  +  ^ 


where  X  is  related  to  the  bulk  diffusivity  of 
gas  in  the  tissue  and  the  blood-tissue  gas 
exchange  half-time. 

Gas  tension  Pd  in  diffusion  region  (=  Pt  in 
tissue  far  way  from  bubble): 


dPd  dPt  Pa-Pt 
dt  dt  t 


3RUT  -  Finite  Tissue  Volume,  Single  Bubble 


Description 

Gas  bubble  is  immediately  surrounded  by 
a  finite,  unstirred,  and  heterogeneously 
perfused  diffusion  region  immersed  in  a 
well-stirred  tissue  compartment.  Diffusion- 
limited  gas  exchange  between  bubble  and 
tissue  occurs  only  in  diffusion  region  of 
fixed  volume. 


Equations 

Rate  of  change  of  bubble  radius: 


dr^ 

dt 


a'tDs 


X  + 


P  —  P  +  ^  Mr 3 

*amb  Mdg  +  ~ 

3r:  3 


Features 

•  Spherically  symmetric  diffusion. 

•  Accounts  for  perfusion  heterogeneity  by 
allowing  spatial  deviations  in  tissue  gas 
tension. 

•  Provides  assessment  of  minimum 
diffusion  region  volume  needed  to 
sustain  bubble  growth. 

Not  applicable  to  multiple  bubbles, 
because  spherical  symmetry  requires  the 
outer  boundary  of  the  diffusion  region  to 
everywhere  adjoin  the  well-stirred 
outermost  tissue  region,  which  does  not 
participate  in  bubble  evolution. 


where  X  is  related  to  the  bulk  diffusivity  of 
gas  in  the  tissue  and  the  blood-tissue  gas 
exchange  half-time. 

Spatial  average  gas  tension  in  the 
diffusion  region: 


dPd 

dt 


1  d 
a'tVd  dt 


(Pi  vo, 


where  Vd  >  Vd(min)  minimum  diffusion 
region  volume. 
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Table  2. ( continued ) 


3RUT-MB  -  Finite  Tissue  Volume,  Multiple  Bubbles 


Description 

Population  of  bubbles  of  different  sizes 
evolving  in  a  finite  volume  of  tissue. 


Equations 

Rate  of  change  of  radius  for  each  bubble 


Features 

•  Bubble  growth  determined  by  extending 
the  3RUT  model  solution  to  a 
spherically  asymmetric  diffusion 
process. 

•  Accommodates  radially  asymmetric 
deviations  in  tissue  gas  tension  caused 
by  perfusion  and  diffusion 
heterogeneity. 

•  Diffusion  region  volume  around  each 
bubble  is  allowed  to  vary 
indeterminately  during  the  bubble 
lifetime. 

•  Accommodates  interactions  between 
multiple  bubbles  by  not  requiring  the 
outer  diffusion  region  volume  of  each 
bubble  to  everywhere  adjoin  the  well- 
stirred  outer  tissue  region. 

•  Prescribes  maximum  bubble  number 
densities  for  a  given  tissue  volume 
based  on  tissue  gas  content. 


dt 


atDs 


A  + 


(Pd-pJ- 


r;  dP. 


amb 


3  dt 


amb 


n  4y  87:  3 

■  Pidg  H - 1 - Mr; 

g  3r;  3 


where  A  is  related  to  X  in  the  2R  and 
3RUT  models  and  to  the  heterogeneity  of 
perfusion  and  gas  diffusivity  in  the  tissue. 

Spatial  average  gas  tension  Pd  in  diffusion 
region  (=  Pt  in  well-stirred  tissue): 


,  N  A 


at  vt 


k=l 


where  N  <  Nmax,  the  maximum  number  of 
bubbles  for  a  given  tissue  volume  Vt. 
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APPENDIX  A.  DERIVATION  OF  MODEL  EQUATIONS 


Assumptions  and  Approximations 

•  Bubbles  are  stationary  and,  therefore,  any  convection  due  to  bubble  movement  may 
be  ignored. 

•  All  gases  involved  in  bubble  evolution  are  ideal. 

•  The  bubble  equilibrates  with  its  surrounding  tissue  much  faster  than  changes  in 
tissue  gas  concentration  due  to  perfusion  effects  or  changes  in  ambient  conditions 
such  as  changes  in  breathing  pressure  and/or  breathing  gas  mixture.  This  allows  the 
diffusion  equation  to  be  solved  ignoring  the  explicit  time-dependent  term  [quasi-static 
approximation.15  The  solution  is  time-dependent  only  through  the  boundary 
conditions  that  it  must  satisfy  (e.g.,  tissue  gas  tension  in  the  well-stirred  region, 
which  may  vary  with  time). 

•  Spherical  symmetry  holds  in  the  diffusion  region.  This  together  with  the  quasi-static 
approximation  leads  to  an  ODE  description  of  the  diffusion  process.  The  assumption 
of  spherical  symmetry  is  relaxed  in  the  final  3RUT-MB  model. 

•  Arterial  blood  is  assumed  to  be  in  equilibrium  with  alveolar  gas,  the  composition  of 
which  is  assumed  to  be  given  by  the  alveolar  gas  equation  13  with  inspired  gas 
assumed  to  consist  of  only  oxygen  (O2)  and  a  single  inert  gas  diluent.  Additionally, 

O2  is  assumed  to  always  be  in  equilibrium  between  tissue  and  bubble,  along  with 
CO2  and  water  vapor.  Resultant  equations  are  readily  elaborated  to  accommodate 
multiple  inert  gas  diluents  and  02  as  a  diffusible  gas  as  described  by  Gerth,  et  al.14 


The  Fick  Relationship  and  General  Expression  for  the  Rate  of  Change  of  Bubble 
Radius 


The  Fick  relationship  is  central  to  the  formulation  of  all  models  of  gas  bubble  dynamics. 
It  relates  the  rate  of  change  in  bubble  gas  content  to  the  gas  concentration  gradient  at 
the  bubble  surface  through  the  gas  diffusivity.  The  Fick  relationship  is  expressed 
quantitatively  as 


1  ci(P ,  V, ) 
RT  dt 


atDsgjAi , 


(Al) 


where  R  is  the  gas  constant,  T  is  temperature,  Pj,  Vi?  and  Aj  are  the  gas  pressure, 
volume,  and  surface  area  of  the  bubble,  respectively,  Ds  and  gi  are  the  gas  diffusivity 
and  gas  tension  gradient  at  the  bubble  surface,  respectively,  and  at  is  the  gas  solubility 
in  tissue  (in  moles  per  unit  volume  per  unit  pressure).  Because  the  product  of  gas 
solubility  and  gas  tension  equals  the  gas  concentration  for  inert  gases,  atgi  is  the  gas 
concentration  gradient  at  the  bubble  surface. 
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Eq.  (A1 )  states  that  the  gas  flux  is  the  product  of  the  permeability  (solubility  x  diffusivity), 
gas  concentration  gradient,  and  surface  area.  For  a  spherically  symmetric  diffusion 
process  around  a  gas  bubble,  which  we  assume  in  all  models  considered  here,  A^Anr? 
and  Vi  =(4n / 3)r;3 ,  where  n  is  the  bubble  radius.  The  bubble  gas  pressure  Pi  is  given  by 


2y  4  7i  3 

P1-Pamb-Pldg+  —  +  —  Mr2 

ri  3 


(A2) 


where  Pamb  is  the  ambient  hydrostatic  pressure,  Pidg  is  the  total  tension  of  gases  that  are 
assumed  to  be  always  in  equilibrium  between  bubble  and  tissue  (e.g.,  metabolic  gases 
and  water  vapor)  and  therefore  to  have  infinite  diffusibilities,  y  is  the  gas-liquid  surface 
tension,  and  M  is  the  tissue  modulus  of  elasticity.3 

Substituting  for  Pi5  Vj,  and  Ai;  we  obtain  the  following  from  Eq.  (A1 ): 


dt 


r 

P 

1  amb 

V 


2v  4tc  3 

_Pidg  + - +  ~Mrf 

r;  3 


An 


(A3) 


where  a't  =atRTis  the  gas  solubility  in  tissue  expressed  in  appropriate  units  including 
the  factor  RT  from  Eq.  (A1 ).  The  left  side  of  Eq.  (A3)  is  expanded  to  obtain 


dP. 


amb 


dt 


- -  +  47TM1-;2 


v  ri 


w 

f  An  3 b 

f 

- V: 

+ 

dt 

/  J 

V  3  'J 

V 

-  P,M„  +  — +  — —  Mr 2  )(47Tr2  =  a'tDsgi  (ati^  ), 

dt 


amb  Mdg 


ri  3 


i.e., 


r,  dP, 


amb 


3  dt 


+ 


f  2  4tt  3 

- +  — Mr.- 

V  3r;  3 


d^ 

dt 


■  + 


p  _  p  2y  An  3  ,  — , _ , 

ramb  Vdg  + - +  ~  Mri 

ri  3 


dr;  ,  _ 

— -atDsgi, 

dt 


i.e., 


p  -  p  ,  8tt  wj 
Vmb  iidg+~  +  ~Mri 


—  =  a(D  g-  amb 
dt  1  S§1  3  dt 


which  is  solved  for  dr/dt  to  yield 


3  The  last  term  in  Eq.  (A2)  is  from  Gernhardt7  and  represents  the  pressure  exerted  on 
the  bubble  by  tissue  elastic  forces,  where  the  elastic  modulus  M  is  given  in  terms  of  the 
bulk  modulus  or  modulus  of  compression,  K,  and  the  mechanically  affected  tissue 
volume,  Vt,  by  M=K/Vt.  As  discussed  in  Appendix  C,  this  term  approximates  expected 
effects  of  tissue  elasticity  on  bubble  pressure  in  some  systems  if  M  is  not  identified  as 
the  elastic  modulus  but  is  instead  treated  as  an  empirical  constant.  Note  that  a  more 
correct  representation  can  readily  be  substituted. 
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^^amb 


(A4) 


dq 

dt 


<Dsgi 


3  dt 


P  —  P  + 

amb  idg 


4y  8n  3 
—  +  — Mr: 
3r  3 


Eq.  (A4)  is  the  general  equation  for  the  rate  of  change  of  bubble  radius  applicable  to  all 
models  of  the  bubble-tissue  system  presented  here.  The  various  models  differ  in  the 
way  this  system  is  represented  to  evaluate  the  pressure  gradient  gj  at  the  bubble 
surface. 


Well-stirred  Tissue  Models 


The  Diffusion  Equation 

In  well-stirred  tissue  models,  gas  exchange  between  bubble  and  a  well-stirred  tissue 
region  occurs  by  diffusion  through  a  diffusion  region  surrounding  the  bubble.  Neglecting 
convection  due  to  bubble  movement,  the  diffusion  process  is  described  without  sources 
or  sinks  by  the  equation 

~=DbV2P ,  (A5) 

at 


where  t  is  time,  P  is  the  diffusion  region  gas  tension,  and  Db  is  the  diffusivity  of  the  gas 
in  the  diffusion  region. 

Assuming  spherical  symmetry  and  denoting  the  radial  distance  from  the  center  of  the 
bubble  by  r,  Eq.  (A5)  becomes 


SP  „  d2P  2  <9P 

— =Dh  — z-+ - 

dt  dr2  r  dr 


(A6) 


We  set  the  time-dependent  term  <3 P/St  to  zero  by  invoking  the  quasi-static 
approximation15  and  reduce  Eq.  (A6)  to 

32P  2  SP 

- r-  + - =  0  . 

dr2  r  dr 


(A7) 


Eq.  (A7)  is  valid  in  the  diffusion  region  n  <  r  <  r0,  where  n  is  the  bubble  radius  and  r0  is 
the  outer  radius  of  the  diffusion  region.  Although  no  longer  explicitly  indicated,  the  gas 
tension  P  still  varies  with  time  because  of  the  time  dependence  of  the  boundary 
conditions  P(n,  t)  =  Pj(t)  and  P(r0,  t)  =  Pt(t),  where  Pi  is  the  bubble  gas  pressure  and  Pt  is 
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the  gas  tension  in  the  well-stirred  region.  The  quasi-static  approximation  allows  solution 
of  Eq.  (A6)  with  separation  of  the  independent  variables  radius  and  time.b 


Solution  of  Diffusion  Equation  and  Expression  for  Pressure  Gradient 

The  general  solution  of  Eq.  (A7)  is  obtained  by  first  solving  for  dP/dr .  Setting  8?/dv  =  z , 
Eq.  (A7)  becomes 

dz  2z 
0r  __T’ 

dz  2 

l.e.,  —  =  — dr, 

z  r 


which  is  integrated  to  yield  ln(z)  =  -ln(r2)  +  K1  ,  where  K-i  is  the  constant  of  integration. 
Equivalently, 


dP_Ki_ 
Sr  ”  r2  ' 


(A8) 


Both  sides  of  Eq.  (A8)  are  integrated  to  yield  P(r)  =  — -  +  K2  .  (A9) 

r 

The  integration  constants  Ki  and  K2  are  evaluated  using  the  boundary  conditions  at  n 
and  r0.  Thus,  K^Pt-PiK^/h  and  K2=(Ptr0-Piri)/h  ,  where  h  =  r0  -  n  is  the  diffusion 
region  thickness.  Substituting  for  Ki  and  K2  in  Eqs.  (A8)  and  (A9)  and  simplifying,  we 
obtain  the  following  expressions  for  the  gas  tension  and  its  gradient  in  the  diffusion 
region,  n  <  r  <  r0: 


P(r,t)  =  Pt 


(pt-p  i)V-r 

h 


(A10) 


and 


b  Eq.  (A6)  can  be  solved  without  the  quasi-static  approximation  by  combining  the  r  and  t 
dependence  of  P  in  a  function  u(r,t)  and  solving  for  P(u). 16,17  However,  the  form 
specified  for  the  function  u(r,t)  determines  how  P  varies  with  t  so  that  such  solutions  are 
not  general.  The  most  recent  attempt  at  such  a  solution17  does  not  account  for  gas  loss 
or  gain  due  to  perfusion  and  is  hence  inapplicable  to  growth  of  a  bubble  in  a  perfused 
tissue  (See  Appendix  B).  A  quasi-static  solution  is  more  practical  as  it  retains  the 
flexibility  to  accommodate  perfusion  heterogeneity,  exercise  effects,  and  multiple 
diffusible  gases. 
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—  =  1  2  '  ,  respectively.  (A1 1 ) 

Note  that  both  the  gas  tension  and  its  gradient  vary  with  time  because  the  boundary 
pressures  Pi  and  Pt  are  functions  of  time. 


I.  Three-Region  Well-Stirred  Tissue  Model  with  Constant  Thickness  Diffusion 
Region  (3RWT-CT) 

Equation  for  Rate  of  Change  of  Bubble  Radius 

From  Eq.  (A1 1 ),  the  pressure  gradient  gi  at  the  bubble  surface  n  is  given  by 


Si 


dP 

dr 


=  (pt-Pi)\  =  (Pt-Pi)^  =  (Pt-Pi) 
r:  h  r:  h 


1  1 

— +  — 
h  riy 


(A12) 


Substituting  the  above  expression  for  gj  in  Eq.  (A4),  we  obtain  the  following  dr/dt 
equation  for  the  well-stirred  tissue  models: 


dt 


ajD; 


f\  l2 

— +  — 

Vh  ri  J 


(Pt-Pi)-T 


^  4Pamb 


3  dt 


P  —  P  + 

L  amb  1  idg  ' 


4y  871  H  3 
—  +  — Mr: 

3r  3 


(A13) 


Tissue  Gas  Tension  in  Well-Stirred  Region 

The  equation  for  determining  the  uniform  gas  tension  in  the  well-stirred  region  is  based 
on  the  expression  for  mass  balance  in  this  region.  Mass  balance  requires  that  the  sum 
of  the  gas  fluxes  and  tissue  gas  contents  in  any  given  time  interval  equal  the  amount  of 
gas  transported  by  perfusion.  Thus,  we  have 

,  roo  r« 

—  JatP(p,t)47rp2dp  =  abQt  J[Pa  -  P(p,t)]47tp2dp  -  f0  ,  (A14) 


where  r^  is  the  outer  radius  of  the  tissue,  otb  is  gas  solubility  in  blood  (in  moles  per  unit 
volume  per  unit  pressure),  Qt  is  blood  flow  per  unit  volume  of  well-stirred  tissue,  Pa  is 
the  arterial  gas  tension,  f0  is  gas  flux  at  the  outer  boundary  of  the  diffusion  region,0  and 
p  is  the  dummy  variable  of  integration  with  respect  to  radial  distance.  Note  that  the  outer 
tissue  radius  is  irrelevant  to  solution  of  either  integral  in  Eq.  (A14),  indicating  that  the 
tissue  need  not  be  spherical  in  shape. 

0  By  convention,  we  take  flux  in  to  be  positive  and  flux  out  to  be  negative. 


A-5 


The  integral  on  the  left  side  of  Eq.  (A14)  evaluates  to  atPt(Vt  -  Vd)  and  the  integral  on  the 
right  side  to  (Pa  -  Pt)(Vt  -  Vd),  where  Vt  is  the  total  volume  of  the  tissue  and  Vd  is  the 
volume  of  the  diffusion  region  given  by 

Vd  =(47t/3)(r03  -ri3)=(47ih/3)(h2  +3rih  +  3ri2).  (A15) 

Substituting  the  evaluated  expressions  for  the  integrals,  dividing  by  atVt  throughout,  and 
rearranging  terms,  Eq.  (A14)  becomes 

£w-v)Pt]  =  (l-v)iz?L- J°  (A16) 

at  x  atVt 

where  v  =  Vd/Vt  is  diffusion  region  fraction  of  the  total  tissue  volume,  and  x=at/abQt  is 

the  blood-tissue  gas  exchange  time  constant.  For  a  large  tissue,  Vt  -»  qo,  v  -»  0,  and  Eq. 
(A1 6)  further  reduces  to  the  familiar  first-order  equation  in  which  the  rate  of  change  of 
tissue  gas  tension  is  proportional  to  the  difference  between  arterial  and  tissue  gas 
tensions. 

Eq.  (A1 6)  requires  evaluation  of  the  gas  flux  f0  at  the  outer  boundary  of  the  diffusion 
region.  The  Fick  relationship  gives  the  gas  flux  f  at  any  radial  distance  r  within  the 
diffusion  region  (n<r<r0)  as  f  =  Dtat4nr2g,  where  g=SP/Sr  is  the  gradient  at  r.  With 
substitution  of  dP/dr  from  Eq.  (All),  r2  cancels  out  and  the  flux  is  seen  to  be 
independent  of  r,  unless  there  is  a  change  in  diffusivity  in  the  diffusion  region,  which  is 
what  we  assume  in  the  3RWT-VTDD  model  discussed  in  a  subsequent  section.  In  the 
3RWT-CT  model,  however,  there  is  no  change  in  diffusivity,  and  the  gas  flux  is  therefore 
the  same  everywhere  in  the  diffusion  region  and  equal  to  the  rate  of  change  in  bubble 
gas  content  given  on  the  left  side  of  Eq.  (A1 ).  Thus, 

f„=f,  =  (A17) 

Ki  dt 

where  fj  is  the  gas  flux  at  the  bubble  surface. 

Using  Eq.  (A17)  to  replace  f0,  and  substituting  a't=atRT  ,  Eq.  (A16)  becomes 


A[(l-V)Pt]=(l-v)^^ 

dt  x 


1  d 

a'tVt  dt 


(Pi  Vi). 


(A18) 


The  time  dependence  of  v  complicates  solution  of  Eq.  (A18).  However,  if  v  is  assumed 
constant  in  each  integration  step,  Eq.  (A18)  becomes 


dl\ 

dt 


P  -P 

ra  rt 


a'tVt[l-v(t)]  dt 


(TO, 


(A19) 
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which  is  the  final  ODE  for  determining  the  tissue  gas  tension  Pt  in  the  well-stirred  region. 
Solution  of  Eq.  (A19)  requires  evaluation  of  v(t)  in  each  integration  step.  On  the  other 
hand,  if  the  diffusion  region  is  negligibly  thin,  the  diffusion  region  volume  Vd  will  be  only 
a  small  fraction  of  the  tissue  volume  Vt.  Thus,  Vd  «  Vt,  and  v  =  Vd/Vt  «  1 .  With  this 
approximation,  the  final  ODE  for  determining  the  tissue  gas  tension  Pt  in  the  well-stirred 
region  is 


dgt  Pa-Pt 

dt  x 


1  d 
a'tVt  dt 


(Pi  Vi). 


(A20) 


The  bubble  radius-time  profile  is  computed  in  the  3RWT-CT  model  using  Eqs.  (A2), 
(A13),  and  (A19)  or  (A20). 

The  radial  gas  flux  in  this  model  is  the  same  everywhere  in  a  thin  diffusion  region  of 
constant  thickness.  As  a  result,  the  diffusion  region  gas  content  remains  unchanged 
while  the  volume  of  the  region  varies  with  bubble  volume.  At  the  same  time,  the  gas 
tension  in  the  region  must  vary  with  the  gas  tension  in  the  well-stirred  region.  Thus  the 
gas  tension  and  gas  concentration  in  the  diffusion  region  vary  in  violation  of  Henry’s 
law,  which  requires  that  the  concentration  of  a  gas  in  any  liquid  region  vary  in  proportion 
to  its  gas  tension. 

It  is  clear  from  Eq.  (A20)  that  bubble  growth  tends  to  deplete  the  gas  tension  Pt, 
reducing  the  driving  force  for,  and  hence  the  rate  of,  bubble  growth.  For  given  t,  Eq. 
(A20)  indicates  that  such  reductions  become  more  severe  as  the  Vt  hosting  a  single 
bubble  decreases. 3  It  is  readily  shown  that  such  reductions  are  equivalent  to  those  that 
occur  as  the  number  of  bubbles  is  appropriately  increased  in  a  given  Vt.  We  begin  by 
rewriting  Eq.  (A14)  for  mass  balance  in  the  well-stirred  region  of  a  tissue  of  volume  Vt 
with  N  bubbles: 


N 


—  J  atP(z,t)dv  =  abQt  J  [Pa  -P(z,t)]dv  -  £(f0)k  , 

d  V*  V*  k=l 


(A21 ) 


where  vts  =  vt  -^(vd)k  is  the  volume  of  the  well-stirred  region  with  (vd)k  being  the 


k=l 

diffusion  region  volume  of  the  kth  bubble,  P(z,t)  is  the  gas  tension  at  any  point  z  in  the 
well-stirred  region,  dv  is  the  elemental  volume  of  integration,  and  (f0)k  is  the  gas  flux  at 

the  outer  boundary  of  the  kth  bubble.  The  integral  on  the  left  side  of  Eq.  (A21 )  evaluates 

f  N  A  (  N 

to  atPt  Vt  ~X(vd)k  ar|d  the  integral  on  the  right  side  to  (pa  -Pt)  Vt  -KU 

V  k=l  )  V  k=l  j 

Substituting  the  evaluated  expressions  for  the  integrals,  dividing  by  otAA  throughout, 
replacing  (f0)k  with  the  rate  of  change  in  gas  content  of  the  kth  bubble,  and  rearranging 

terms,  Eq.  (A21)  becomes 
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1  N  A 

1  d 


^[(1-v)P,]  =  (1-v)5lA 

dt  T  atVt  dt 


(A22) 


where  v  is  now  k=1 


vt 


.  If  v  remains  negligibly  small,  Eq.  (A22)  reduces  to  the  multiple 


bubble  analog  of  Eq.  (A20): 

dt  i  a;vt^dtliJk 


(A23) 


It  follows  from  Eq.  (A20)  and  (A23)  that  a  system  with  N  identical  bubbles 

f  N  d  j  A 

^  —  (p; Vi  )k  =  N  —  (PjVj )  in  tissue  volume  (vt)N>1  will  be  identical  to  a  system  with  only 

v  k=1  ^  ^t 

a  single  bubble  in  a  tissue  volume  (vt)N=1  given  by 

(v,)„ 


tV,)N.,  = 


/N>1 


N 


(A24) 


Eq.  (A24)  is  inverted  to  show  that  equivalence  of  the  two  systems  will  be  maintained 

1  N 

provided  they  have  equivalent  bubble  number  densities:  - — r —  =  7 — - — . 

(v,L„  (v,)N>, 


Numerical  Procedure  for  Solving  Eq.  (A20) 
Let  X  =  PiVi/RT.  Then  Eq.  (A20)  becomes 


dPl+PL=Pa  l  dX 


dt  x  x  atVt  dt 


(A25) 


Both  sides  of  Eq.  (A25)  are  multiplied  by  the  integrating  factor  et/x  to  obtain 

1  dX 
atVt  dt 


f  (40  o  A 

dt  x 


|  gt/x  _  Pa  ct/x  ^  dX  ct/T 


y 


i.e.,  —  (pte^T)=^-e^ - - — . 

dt v  x  atVt  dt 


(A26) 
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Letting  Pa(s)  =  Pa(t)  +  Pa(s-t)  in  the  interval  t  <  s  <  t  +  At  with  fixed  values  for  Pa(t)  and  Pa , 

and  integrating  between  t  and  t  +  At,  where  At  is  the  integration  step  size  and  s  is  the 
dummy  variable  of  integration,  we  obtain  the  following  from  Eq.  (A26): 

Pt  (t  +  At)e(t+At)/x  -  Pt  (t)et/x 


t+At 


t+At 


dX 


-  f[Pa(t)  +  Pa(s-t)^ds-—  f  es/x  — 
t  Jt  atVt  Jt  dt 


ds 


Pa(t) 


t+At  t+At 


t+At 


,s/x 


dX 


L  I  l  y.  I  I  ..  l  I  !— 

Law  f eS//xds  +  —  f (s - t)es/Tds - —  f  e  — 

J  t  J  a  V,  J  dt 

t  t  ^t  vt  t 


Pa  (t) 


ds 


g(t+At)/x  _  e t / x 

+  Pa 

(s-t)es/x-xes/\ 

,  t+At 

t+At  1  ,  gS/TdX 

1  a'tVt  ]  dt 

ds 


t+At 


s/x 


dX 


Pa(t)[e(t+At)/T-et/T]  +  Pa[(At)e(t+At)/T-Te(t+At)/x+Tet/xl - —  f  e 

cxtVt  Jt  dt 


ds.  (A27) 


Multiplying  both  sides  of  Eq.  (A27)  by  e  (t+AtVT; 

Pt (t  +  At) - Pt (t)e_At 7 x  =  Pa[l-e_At/x]  +  Pa(At)-Pax[l-eAt/xl  6 

and  therefore, 


~(t+At)/x 


t+At 


s/x 


atVt  Jt 


dX 

dt 


ds , 


Pt  (t  +  At)=Pt  (t)e_At/x  +  Pa  [1  -e“At/x]  +  Pa  (At)  -  Pax|l-e 


3  At/x 


-(t+At)/x 


t+At 


s/x 


atYt  Jt 


dX 

dt 


ds 


Pt  (t)  -Pt  (t)[l  -e-At/x]  +  Pa  [l-e_At/x]  +  Pa  (At)  -  Pax|l-e 


sAt/x 


-(t+At)/x 

atVt 


t+At 


s/x 


dX 

dt 


ds 


=Pt(t)  +  Pa(At)  +  (pa -Pax-Pt)tl-e 


-At/xn  6 


-(t+At)/x 


t+At 


s/x 


atV,  Jt 


dX 

dt 


ds. 


(A28) 
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The  integral  in  Eq.  (A28)  can  be  evaluated  by  approximating  X(s)  by  a  straight  line  in 

the  interval  [t,  t  +  At],  i.e.,  by  letting  X(s)=X(t)+X(t  +  At^  X(t\s-t)  for  t  <  s  <  t+At.  Then,  — 

At  dt 

the  integrand  of  the  integral  term  on  the  right  side  of  Eq.  (A28)  is  simply  the  slope 

X(t  +  At)-X(t) 


in 


At 


,  and 


-(t+At)/x 


t+At 


,s/x 


atVt  Jt 


dX 

dt 


ds=e-^"/'X(,  +  A,)-X(,)'J;„,ds 


atVt 


At 


t+At 


r(t+At)/x  X(t  +  At)-X(t) 
atVt  At  J 


s-(t+At)/x 


X(t  +  At)-X(t) 


atVt 


At 


e(t+At)/x  _ e t/x 


X(t  + At)-X(t)  T 


atVt  At 


l-e_At/T 


(A29) 


Substitution  of  the  expression  for  the  integral  given  by  Eq.  (A29)  into  Eq.  (A28)  yields 


Pt  (t  +  At)  =  Pt  (t)  +  Pa  (At)  +  (pa -Pax- Pt  \\  -  e“At/T  ] 


-At  /  X  n  X(t  + At)-X(t)  X 


cqVt  At 


l-e“At/x 


=  Pt(t)  +  Pa(At)  + 


Pa  -Pa^-Pt  - 


X(t  +  At)  -  X(t)  X 


atVt 


At 


1  —  e_At/T 


(A30) 


With  a  steady  arterial  tension,  pa  =  o,  and  Eq.  (A30)  reduces  to  the  expression  for 
approximating  Pt  given  in  reference  3.  Then  as  Vt  -»  qo,  X(t  +  At)-X(t)  anc|  (A30) 


cxtVt 


further  reduces  to 

Pt(t  + At)=Pt(t)  +  (Pa  —  Pt  )[l  —  e_At^T  ] , 


i.e.,  Pt (t  +  At)- Pt (t)  =  (Pa  -Pt)[l- 


i.e., 


Pt  (t  +  At)— Pt  (t) 
At 


=  (Pa  —  Pt  ) 


e 

1-e 


-At/x 


-At/x 


At 
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as  At  — >  0, 


(A31 ) 


which  is  the  same  as  Eq.  (A20)  without  the  second  term  on  the  right  side  that  accounts 
for  gas  loss  into  the  bubble  or  gain  from  the  bubble.  In  applying  Eq.  (A30),  we  first  solve 
Eq.  (A13)  for  n  assuming  Pt  to  be  constant  in  the  interval  [t,  t  +  At],  in  accord  with  the 
quasi-static  approximation.  We  then  use  Eq.  (A30)  to  update  Pt.  Thus,  Eq.  (A30)  needs 
to  be  used  only  once  during  each  integration  step.  Also,  other  approximations  for  X(s)  in 
the  interval  [t,  t  +  At],  e.g.,  exponential,  can  be  used  to  obtain  expressions  similar  to  Eq. 
(A30)  for  Pt(t+At). 


II.  Three-Region  Well-stirred  Tissue  Model  with  Varying  Thickness  Diffusion 
Region  (3RWT-VTDD) 

The  violation  of  Henry’s  law  that  occurs  in  the  diffusion  region  of  a  growing  or  resolving 
bubble  in  the  3RWT  model  can  be  avoided  if  the  diffusion  region  thickness,  and  hence 
its  volume,  is  varied  so  as  to  produce  changes  in  diffusion  region  gas  concentration  and 
tension  appropriate  for  maintaining  a  constant  gas  content. 


Equation  for  Rate  of  Change  of  Bubble  Radius 

The  equation  for  the  rate  of  change  of  bubble  radius  is  the  same  in  this  case  as  in  the 
3RWT-CT  model.  Eq.  (A13)  for  dr/dt  holds  with  the  difference  that  the  diffusion  region 
thickness  h  varies  with  time  and  must  be  computed  at  each  integration  step.  The 
equations  required  for  these  updates  of  h  are  derived  below. 


Equations  for  Diffusion  Region  Thickness 

ro 

The  gas  content  of  the  diffusion  region  is  Ud=JatP(p,t)47tp2dp .  (A32) 

ri 

Evaluating  the  integral  in  Eq.  (A32)  using  P(p,t)  from  Eq.  (A10),  we  obtain 


Ud=<xtjPt  47tp2dp  -  at  4n(Pt  -  P; )  j-  j  (rDp  -  p2  )dp 


atPtVd-at47r(Pt  -  Pj) 7- 

h 


r  (r2  —  r-2 )  r3  —  r3 

oV  o  1  /  o  1 
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:atPtVd-at^(Pt-Pi)ri(r°  fi) 
3  h 


3r0(ro  +fi)-2(r02  +r0r;  +rf) 


2tc  o 

a  t  Pt ' Vd  -  a  t  y  (Pt  -  Pj  )r,  [r0  (r0  +  r; )  -  2r, 


?  71  -> 

atptVd~at  y(pt  -pi)ri[(ri  +  h)(2ri  +  h)  -  2rf 


■  a. 


Pt^d  +h)(Pt  -  Pi) 


(A3  3) 


Substituting  for  Vd  from  Eq.  (A15)  and  for  Pj  from  Eq.  (A2),  Eq.  (A33)  becomes 


TT  2lZ 

ud=yat 


2Pt  (r03 -r,3 )  +  rih(3ri  +h)f P^mb  3+^ Mr3  -Pt 


v  ri  3 


(A34) 


where  p;mb  =Pamb  -Pidg.  Substituting  r0  =  n  +  h  and  rearranging  terms  in  Eq.  (A34),  we 
obtain 


2Pt|(r1  +  h)3-r3|+r1h(3r1+h)! 


Pamb  +  ~  +  ~  Mr  3  -  Pt 
ri  3 


2 n  at 


i  e.,  (2Pt)[3rrh  +  3r1h2  +  h3j+(3r1h  +  hi)^P'mbr1  +  2y +— Mr,4  -  Ptr, 


4tx 


:_3_Ud 

2 n  at 


2  ,  A  u 3  _  3  Ud 


l.e.,  A1h+A,h-+A3hJ= - -, 

2n  at 


(A3  5) 


where,  A{  =  6yrj+3(Pt  +  p'mb)rf  +4nMri5 , 

Ajr  . 

A2=r1(5Pt  +  P;mb)  +  2y  +  — M^4,and 
A3  =  2Pt. 

We  now  argue  that  the  gas  content  of  the  diffusion  region  as  given  by  Eq.  (A35)  can 
neither  be  zero  nor  a  constant.  If  Ud  =  0,  then  h  =  0  is  the  only  admissible  solution  of  Eq. 
(A35).  Any  other  admissible  value  of  h  would  have  to  be  a  positive  real  root  of  the 
quadratic  equation,  A1+A2h+A3h2=0 ,  which  is  not  possible;  the  two  possible  real  roots 
are  negative  because  the  coefficients  A-i,  A2,  and  A3  are  all  greater  than  zero.  If  Ud  is  a 
constant,  the  gases  in  the  diffusion  region  follow  Henry’s  law,  but  the  thickness  h  rapidly 
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diminishes  as  the  bubble  grows  (Fig.  2,  Ref.  4).  With  constant  gas  diffusivity,  the 
diffusion  region  offers  decreasing  diffusive  resistance  to  gas  flux  into  the  bubble.  As  a 
result,  gas  tensions  in  the  bubble  and  its  surroundings  rapidly  approach  equilibrium,  and 
the  dynamics  of  continued  bubble  growth  are  no  longer  diffusion-limited.  Once  the 
diffusion  region  is  brought  into  theoretical  conformance  with  Henry’s  law,  bubble  growth 
can  evidently  remain  diffusion  limited  only  if  the  gas  content  of  the  diffusion  region  is 
allowed  to  vary.  Given  that  the  diffusion  region  is  not  perfused,  such  variation  can  occur 
only  if  the  fluxes,  f0  and  fh  are  unequal. 

Requisite  flux  inequality  is  obtained  if  the  diffusivity  at  the  bubble  surface  Ds  is  different 
from  the  diffusivity  Db  in  the  remainder  of  the  diffusion  region.  Accordingly,  we  assume  a 
higher  resistance  to  diffusion  at  an  infinitesimally  thin  bubble  surface  with  diffusivity  Ds 
less  than  Db.  Because  the  concentration  gradient  and  surface  area  at  such  a  surface 
are  the  same  for  both  f0  and  fj,  we  find  by  applying  the  Fick  relationship  that 

Il=  Dsg.gA  =£l.  (A36) 

fo  DbatgoAo  Db 


Thus,  a  differential  diffusivity  at  the  bubble  surface  with  constant  values  of  Ds  and  Db 
results  in  proportionality  of  fluxes  fj  and  f0.  Using  Eq.  (A36)  and  expressing  the  flux  fj  as 
the  rate  of  change  of  bubble  gas  content,  we  obtain  the  following  mass  balance 
equation  for  the  diffusion  region: 


dUd 

dt 


=f0-fi  =(P-1) 


1 

RT 


d(Pj  Vj ) 


(A37) 


where  p  =  Db/Ds  is  the  diffusivity  ratio  (p  >  1  since  Ds  <  Db).  Eq.  (A37)  is  integrated  to 
yield 

Ud=atkdPiVi ,  (A3  8) 

where  kd=(p-l)/a't .  Eq.  (A37)  permits  a  zero  value  for  the  constant  of  integration  so 
long  as  p  >  1 .  Accordingly,  we  have  taken  the  constant  of  integration  to  be  zero  so  that 
h  vanishes  along  with  n  upon  bubble  resolution.  In  the  degenerate  case  of  p  =  1 ,  Eq. 
(A37)  reduces  to  dUd/dt  =  0,  and  Ud  is  the  non-zero  constant  content  characteristic  of 
the  3RWT-CT  model.  This  is  consistent  with  Eq.  (A35)  in  which  Ud  cannot  be  zero. 


We  obtain  the  following  cubic  equation  in  h  by  substituting  Ud  from  Eq.  (A38)  in  Eq. 
(A35): 


A1h+A2h2+A3h3  =(3 / 27i)kdP1Vi .  (A39) 

Recalling  the  definitions  of  the  coefficients  A-i,  A2,  and  A3  given  with  Eq.  (A35),  Eq. 
(A39)  indicates  that  the  diffusion  region  thickness  is  a  function  of  the  differential 
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diffusivity  ratio,  the  size  and  content  of  the  gas  bubble,  and  the  gas  tension  in  the  well- 
stirred  tissue  region.  The  diffusion  region  volume  is  then  given  by  Eq.  (A15). 


Tissue  Gas  Tension  in  Well-Stirred  Region 

Except  for  the  fact  that  f0  =  pfi  with  p^1 ,  the  equation  for  the  uniform  gas  tension  Pt  of 
the  well-stirred  region  is  the  same  here  as  in  the  3RWT-CT  model.  Accordingly,  Eq. 
(A17)  is  modified  as  follows: 

fo=Pfi  =  ^PiVi)-  (A40) 

RT  dt 

Eq.  (A40)  leads  to  the  following  final  equation  for  determining  Pt,  an  equation  that 
corresponds  to  Eq.  (A18)for  the  3RWT-CT  model: 

|[<i  -  v)p,i = <a41  ) 

dt  x  atVt  dt 

The  diffusion  region  volume  fraction  v  cannot  be  considered  negligibly  small  in  the 
3RWT-VTDD  model.  The  solution  of  Eq.  (A41)  is  involved  because  of  the  presence  of 
the  dv/dt  term  on  the  left  side.  However,  Eq.  (A41)  can  be  simplified  to  compute  Pt  by 
assuming  v  to  be  a  constant  in  each  integration  step.  If  v  is  a  constant,  Eq.  (A41 ) 
simplifies  to 

=  - P - A(pVj).  (A42) 

dt  t  a'tVt[l -v(t)]  dt 


Eq.  (A41)  is  modified  for  a  system  with  N  bubbles  as  follows: 


>-v)Pt]  =  (l-v)^ 
dt 


-Pt 


(A43) 


where  v  = 


I(vd)k 


k=l 


Vt 


with  (vd)k  evaluated  for  each  bubble  using  Eqs.  (A39)  and  (A15). 


As  above,  v  cannot  be  neglected,  but  if  it  is  again  assumed  to  be  constant  in  each 
integration  step,  Eq.  (A43)  reduces  to  the  multiple  bubble  analog  of  Eq.  (A42): 


dPt 

dt 


P 


t  a;vt(i-v(t))sdt 


k  ■ 


(A44) 
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Eqs.  (A2),  (A13),  (A39)  and  (A42)  [or  Eqs.  (A2),  (A13),  (A39),  one  for  each  bubble,  and 
Eq.  (A44)]  are  the  requisite  equations  for  computing  the  bubble  radius-time  profile  in  the 
3RWT-VTDD  model.  While  Eqs.  (A2)  and  (A13)  are  the  same  as  for  the  3RWT-CT 
model,  Eq.  (A42)  replaces  Eq.  (A20)  and  Eq.  (A44)  replaces  Eq.  (A23)  for  computing 
the  gas  tension  in  the  well-stirred  region.  Solutions  of  Eq.  (A42)  or  Eq.  (A44)  are  more 
involved  than  the  solutions  of  the  corresponding  Eqs.  (A20)  or  (A23)  as  they  remain 
valid  in  only  small  intervals  of  t.  The  numerical  solution  given  by  Eq.  (A30)  for  the 
3RWT-CT  model  is  applicable  to  Eq.  (A42)  if  v(t)  remains  at  a  steady  value  in  each 
integration  interval  and  p/vt [l— v(t)]  is  substituted  for  i/vt  .  To  avoid  large  changes  in  v(t) 
in  any  integration  interval,  it  may  be  necessary  to  use  a  smaller  step  size  than  needed 
for  the  3RWT-CT  model.  Also,  the  3RWT-VTDD  model  requires  additional  computations 
to  solve  Eq.  (A39)  for  the  varying  diffusion  region  thickness  around  each  bubble.  Note 
that  when  the  model  is  elaborated  to  accommodate  multiple  diffusible  gases,  the 
diffusion  region  thickness  and  volume  is  different  for  each  gas  due  to  the  different 
diffusivity  ratios,  tissue  gas  tensions,  and  bubble  pressures  of  each  gas. 


Unstirred  Tissue  Models 


The  Diffusion  Equation 

All  unstirred  tissue  models  are  based  on  solution  of  the  spherically  symmetric  diffusion 
equation  given  by  Eq.  (A6),  but  with  an  added  sink  term  to  accommodate  the  effects  of 
perfusion  on  bubble  evolution.  Accordingly,  Eq.  (A6)  takes  the  following  form  for  these 
models: 


SP  n 

at  TT  =  atDb 

fit 


S2P  +  2  0P 
fir2  r  fir 


~abQt(p-ps)> 


(A45) 


where  Qt  is  blood  flow  per  unit  tissue  volume,  at  and  ab  are  gas  solubilities  in  tissue  and 
blood,  respectively  (in  moles  per  unit  volume  per  unit  pressure),  and  Ps  is  the  sink 
pressure. 

Note  that  all  quantities  in  Eq.  (A45)  are  expressed  in  terms  of  gas  content  rather  than 
gas  tension.  The  pressure  difference  (P  -  Ps )  acts  as  a  sink  if  P  >  Ps  and  as  a  source  if 
P  <  Ps.  Invoking  the  quasi-static  approximation,  we  set  dP/dt  to  zero,  and  reduce  Eq. 
(A45)  to 


atDb 


fi2P  2  5P 

— T+“T- 

fir2  r  fir 


-abQt(p_ps)  - 


i.e., 


fi2P  2  0P  2 

— -+ - X 

fir  r  fir 


(p-ps)=0, 


(A46) 
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where  ?.2=abQt/atDb . 

Strictly  speaking,  the  sink  pressure  Ps  equals  the  arterial  gas  tension  Pa  at  locations  in 
the  diffusion  region  where  capillaries  are  present  and  is  zero  at  all  other  locations.  Also, 
the  blood  flow  may  not  be  the  same  in  all  capillaries.  For  mathematical  tractability,  we 
assume  a  uniform  blood  flow  Qt  in  the  diffusion  region  and  a  sink  pressure  Ps  that  is 
different  from  Pa.  The  sink  pressure  is  a  constant  independent  of  position  in  the  2R 
model,  is  a  function  of  only  the  radial  distance  in  the  3RUT  model,  and  is  position- 
dependent  in  the  3RUT-MB  model. 


III.  Two-Region  Model  (2R) 

Solution  of  Diffusion  Equation  and  Expression  for  Pressure  Gradient 

Eq.  (A46)  is  solved  in  the  2R  model  using  the  following  boundary  conditions:  P(n,  t)  = 
Pi(t),  the  bubble  gas  pressure,  and  the  gradient  dP/dr  vanishes  at  large  distances  from 
the  bubble.  The  latter  condition  ensures  that  tissue  regions  far  away  from  the  bubble  do 
not  contribute  to  bubble  evolution. 


Using  the  fact  that  Ps  is  independent  of  r,  we  re-write  Eq.  (A46)  as 


d2V  25P' 
dr2  r  dr 


-12P'  =  0, 


(A47) 


where  P'=P-PS. 


Define  P'=- .  Then, 
r 

5P'  l  dz  z 
dr  r  dr  r2 


and 


a2p' _ia2z  2  dz  1 2z_ia2z  2ap' 
ar2  r  ar2  r2  Sr  r3  r  dr2  r  dr 

In  terms  of  the  variable  z,  Eq.  (A47)  becomes 

a2  z  o 

— -^--A.2z  =  0  .  (A48) 

dr2 
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The  parameter  X  >  0,  and  therefore,  the  general  solution  to  Eq.  (A48)  is: 
z(r)=K1e~>Lr+K2e+^r ,  where  K-i  and  K2  are  integration  constants.  Replacing  z  by  P'r  and 
then  P'  by  (P-Ps),  we  obtain 

P(r)=Ps  +— e_x,r+— ^-e+A,r  _  (A49) 

r  r 

The  constant  K2  associated  with  the  growing  exponential  term  is  zero  because  of  the 
condition  that  the  gradient  dP/dr  -»  0  as  r  -»  qo.  The  constant  K-i  is  evaluated  using  the 

boundary  condition  at  n.  We  have,  P^Pj+^-e-^  .  Therefore,  K1=(Pi-Ps)rie+Aj‘ ,  and  the 

ri 

final  equation  for  P(r,  t)  is 


P(r,t)=Ps-(P8-Pj)-!-e 

r 


P  —  A.(r— r; ) 


(A50) 


The  second  term  on  the  right  side  of  Eq.  (A50)  vanishes  and  P(r,  t)  equals  Ps  as  r  -»  qo. 
This  means  that  the  sink  pressure  Ps  in  the  2R  model  is  the  tissue  gas  tension  Pt  far 
away  from  the  bubble.  The  time-dependency  of  P(r,  t)  arises  from  changes  in  Pt  with 
time  brought  about  by  blood  perfusion. 

Setting  Ps  =  Pt  in  Eq.  (A50),  we  obtain  the  following  expression  for  the  pressure  gradient 
by  direct  differentiation: 


<9P 

dr 


-=(pt-pi) 


rfx+lL 

r  r2 


(A51 ) 


Equation  for  Rate  of  Change  of  Bubble  Radius 

From  Eq.  (A51 ),  the  pressure  gradient  at  the  bubble  surface  gj  is  given  by: 


gi  = 


SP 

dr 


=(pt-pi) 


X-\ — 


V 


xi  J 


(A52) 


Substituting  the  above  expression  for  gi  in  Eq.  (A4),  we  obtain  the  following  dn/dt 
equation  for  the  2R  model: 


dr, 

dt 


a'tDs 


X  + 


(Pt  —  Pi ) 


r,  dP. 


amb 


"i  J 


dt 


P  —  P  ^  |  Mr3 

^amb  ^idg  ^  3j-.  ^  3  '  ^ 


(A53) 
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It  is  not  necessary  that  the  bubble  surface  diffusivity  Ds  be  different  from  the  tissue  bulk 
diffusivity  Db  in  this  model.  This  is  because  the  radius-time  profile  corresponding  to 
different  values  of  Ds  and  Db  can  be  obtained  by  setting  Ds  =  Db  and  adjusting  other 
parameters  of  the  model  without  loss  of  generality. 


Tissue  Gas  Tension  in  Far  Region 

Eq.  (A51 )  indicates  that  dP/dr  is  zero  only  when  r  reaches  infinity,  from  which  it  follows 
that  the  tissue  volume  involved  in  bubble-tissue  gas  exchange  is  theoretically  infinite  in 
the  2R  model.  As  the  tissue  volume  Vt  approaches  infinity,  the  term  containing  Vt 
vanishes  in  Eqs.  (A20)  and  (A42)  for  determining  the  tissue  gas  tension  Pt.  The  tension 
far  away  from  the  bubble  is  unaffected  by  the  gas  losses  and  gains  by  the  bubble,  and 
is  given  by: 

dP  P  -P 

.  (A54) 

dt  x 

Eqs.  (A2),  (A53),  and  (A54)  constitute  the  2R  model,  applicable  for  determining  the 
growth  of  a  single  bubble  immersed  in  a  tissue  of  infinite  size.  They  involve  the  tissue 
gas  tension  Pt  far  away  from  the  bubble.  Although  one  would  intuitively  discount  any 
role  for  this  gas  tension  in  determining  bubble  growth,  it  turns  out  that  Pt  is  the  same  as 
the  spatial  average  tissue  gas  tension  in  the  infinite-sized  diffusion  region.  This  is  shown 
by  considering  the  spatial  average  tissue  gas  tension  Pd  between  r  =  n  and  r  =  r0  in  the 
diffusion  region  having  a  volume  Vd.  From  Eq.  (A50)  with  the  sink  pressure  Ps  =Pt ,  Pd  is 
obtained  as 


Ao  Ao  Ao 

| P(p)4np2dp  =  | Pt47tp2dp  —  J* (Pt  —  Pj  )e~>(p  r| Tj  4jipdp  , 

ri  ri  ri 


i.e. 


VdPd  = 


VdPt-j(Pt-Pi)e- 


Mp-r;) 


r;  4npdp 


(A55) 


The  second  term  on  the  right  side  of  Eq.  (A55)  is  finite,  and  therefore,  Pd  -»Pt  as 

Vd  ->  oo.  The  sink  pressure  Pt  determined  by  Eq.  (A54)  describing  the  perfusion  process 
is  consistent  with  the  diffusion  equation  for  this  model.  To  show  this,  we  first  modify  Eq. 


(A45)  by  dividing  both  sides  by  at  and  substituting 


S2P  +  2  8? 
dr2  r  dr 


J__S  2  5P 
r2  fir  _  fir 


and 


ab^‘  =-  to  obtain 

at  t 
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(A56) 


3P  13  2  3P 

—  =  Dh  — —  r  — 
3t  r2  5r  [  3r 


P-P 

x  S 


The  spatial  average  tissue  gas  tension  Pd  between  r  =  n  and  r  =  r0  in  a  diffusion  region 
of  finite  volume  Vd  is  obtained  by  integrating  Eq.  (A56)  with  respect  to  the  radial 
distance.  Thus, 


—  jp(p)47tp2dp  =  Dt 


.  7  3P 

4nr  — 
3r 


1  r° 

-J[P(p)-Ps]47ip2dp, 


I.e.,  Vd 


dPd  f0-fi  „  Pd-Ps 


dt  at 


V, 


i.e., 


dPd  f0-fi  Pd-Ps 


dt  atVd  x 


(A57) 


where  f0  and  are  the  gas  fluxes  at  r  =  r0  and  r  =  n  respectively.  As  Vd  ->  qo,  the  first 
term  on  the  right  side  of  Eq.  (A57)  drops  out,  and  the  spatial  average  tissue  gas  tension 
Pd  ->Pt.  Then,  Eq.  (A57)  reduces  to  Eq.  (A54)  with  Ps  =  Pa.  Thus,  the  sink  pressure  is 
the  arterial  pressure  Pa  in  the  diffusion  equation  without  the  quasi-static  approximation, 
i.e.,  in  the  diffusion  equation  with  the  time-dependent  term  3P/3t .  However,  it  is  hard  to 
interpret  Eq.  (A56)  with  Ps  =  Pa,  which  implies  presence  of  capillary  blood  at  every  point 
within  the  diffusion  region. 


With  the  quasi-static  approximation,  3P/3t  is  set  to  zero  and  the  dPd/dt  term  on  the  left 
side  of  Eq.  (A57)  is  absent.  Then,  as  Vd  ->  oo  ,  pd  ->pt,  and  Ps  =  Pt-  In  this  case,  the  sink 
pressure  is  the  uniform  gas  tension  Pt  prevailing  in  the  infinitely  large  tissue  volume. 


IV.  Three-Region  Unstirred  Tissue  Model  (3RUT) 

In  this  model,  the  sink  pressure  varies  with  radial  distance  r  according  to  the  equation, 

Ps(r)  =  Pd  +  F(r),  (A58) 

where  Pd  is  the  uniform  component  of  Ps(r)  independent  of  r  and  equal  to  the  spatial 
average  gas  tension  in  the  diffusion  region, d  and  F(r)  is  any  function  of  r.  The  function 
F(r)  is  the  non-uniform  component  of  Ps(r)  that  represents  radial  deviations  in  sink 
pressure  from  its  uniform  component.  We  assume  these  deviations  to  be  the  same  at  all 
points  lying  on  a  sphere  of  radius  r  from  the  center  of  the  bubble,  consistent  with  our 

d  The  "bar"  notation  to  indicate  spatial  average  was  omitted  in  the  original  description  of 
this  model  in  reference  5. 
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assumption  of  spherical  symmetry  for  describing  the  diffusion  process.  It  should  be 
noted  that  F(r),  and  hence  Ps(r),  are  functions  of  time  as  well.  However,  to  simplify  the 
notation  we  will  not  explicitly  indicate  the  time-dependence  of  F(r)  and  its  derivatives 
and  integrals. 


Solution  of  Diffusion  Equation 

With  Ps  represented  by  Eq.  (A58),  the  diffusion  Eq.  (A45)  becomes 


5P  ^ 

“t  ^r  =  atDb 

fit 


fi2P 

dr2 


2fiP 
r  fir 


ab 


Qt[p-Pd-F(r)], 


which,  under  the  quasi-static  approximation,  reduces  to 


fi2P  |  2d? 
dr 2  r  dr 


-  A,2[p  -  Pd  -  F(r)]=  0  , 


(A59) 


(A60) 


where  i2  =abQt/atDb,  as  defined  for  Eq.  (A46). 

The  transformation  P  =  u/r  yields 

5P  1  fiu  u  ,  fi2P  1  82 u  2  fiu  2u 

—  = - ,  and  — -  = - - - - —  +  — 

Sr  r  fir  r2  fir2  r  dr2  r2  fir  r3 


and  transforms  Eq.  (A60)  to 


1  S2u  2  fiu  2u 

2 

1  fiu  u 

-  ? 

u  —  _ .  . 

- - _|_  - 

+  — 

— 

-Id 

- Pd  -  F(r) 

r  fir2  r2  5r  r3 

r 

_r  fir  r2  _ 

_  r 

=  0, 


ie.. 


r  fir2 


--Pd-F(r) 

r 


=0, 


i.e., 


fi2u 

fir2 


-^2u  =  -A.2r[pd  +F(r)]. 


(A61 ) 


Eq.  (A61 )  is  a  linear  differential  equation  involving  a  single  independent  variable  r.  The 
general  solution  of  Eq.  (A61)  is 

u(r)=Kie“Xr+K2eXr+v(r),  (A62) 

where  the  first  two  terms  with  associated  constants  Ki  and  K2  constitute  the 
homogeneous  part,  and  \|/(r)  is  the  particular  integral  resulting  from  the  function  on  the 
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right  side  of  Eq.  (A61 ).  We  determine  the  particular  integral  by  the  method  of  variation  of 
parameters.18  In  this  method,  we  write  \\>(r)  as  a  combination  of  the  two  independent 
solutions  that  form  the  homogeneous  part.  That  is, 

V(r)=Ui(r)e^r+u2(r)eXr  .  (A63) 


We  determine  the  functions  ui(r)  and  u2(r)  to  satisfy  Eq.  (A61 )  with  the  requirement  that 
the  second  derivatives  of  Ui  and  u2  be  absent.  This  requirement  yields  first-order 
differential  equations  for  the  unknown  functions  u-i(r)  and  u2(r)  and  leads  to  their  integral 
solutions  as  indicated  below,  where  we  ignore  time-dependence  and  use  the  notation 
for  total  derivatives  in  place  of  partial  derivatives. 


Differentiating  \j/(r)  with  respect  to  r,  we  obtain  from  Eq.  (A63) 


dv|/ 

d7 


(-A)Ule^r  +(A)u2eXr 

+ 

dui  -Xr  du  2  Xr 
- e  h - e 

dr  dr 

(A64) 


In  order  that  second  derivatives  of  ui  and  u2  be  absent  when  v|/(r)  is  used  to  satsisfy  Eq. 
(A61 ),  we  equate  the  second  term  in  brackets  on  the  right  side  of  Eq.  (A64)  to  zero, 


du,  e-i, +  dulel,=0  _ 


dr 


dr 


(A65) 


and  obtain 


d\| / 

d7 


(-A)Ule“Xr  +(X)u2eXr 


and 


d\ 

d^ 


,2  -Xr  0  4U|  _^r 

A  Uie  -A - e 

dr 


1 2  Xr  ,  du2  Xr 

A  u2e  +  A - e 

dr 


(A66) 


Since  \\ /(r)  must  satisfy  Eq.  (A61 ),  we  use  Eqs.  (A63)  and  (A66)  in  Eq.  (A61 )  to  obtain 


1 2  -Xr  dUj  _Xr 

A  u  i  e  —A - e 

1  dr 


^  2  Xr  .  ^  2  _  Xr 

A  u2e  +A - e 

2  dr 


-X2 


u1e_A'r+u2eXr 


-A2[pdr+F(r)], 


du. 


du- 


i  e.,  — ^-e  >r  -^Te>r  =  Ar[pd  +F(r)]. 

dr  dr 


(A67) 


Eqs.  (A65)  and  (A67)  are  two  simultaneous  equations  in  du-i/dr  and  du2/dr,  from  which 
we  obtain 


du 


*  ^+F(r)]2 


A,r 
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and 


du2 

dr 


■  =  -Aj[pd+F(r)]- 


-Xr 


Xp 


which  yield  uj(r)=  fAp[pd  +F(p)]- — dp 


(A68) 


and 


Xp 


u2(r)  =  -J  7p[pd+F(p)]^— dp. 


(A69) 


Substituting  the  above  expressions  for  ui(r)  and  U2(r)  in  Eq.  (A65),  we  obtain 

¥(r)  =  e-u|rA.p[pd  +F(p)]^dp-eXr|r^p[pd  +F(p)]^dp 


J^p[pd+F(p)] 


eX(p-r)_e-X(p-r) 


dp 


-  X?. 


'd  J  p  sinh{A,(p  -  r)}dp  +  X J  pF(p)  sinh{A,(p  -  r)}dp  . 


(A70) 


Integrating  the  first  of  the  two  integrals  in  Eq.  (A70)  by  parts  yields 


v(r)  =  7Pd 


pcosh{^(p  -r)}  sinh{A,(p -r)} 


X 


X1 


p=r  r 

+  x|pF(p)sinh{X,(p -r)}dp 


=  Pdr  +  7G(r) , 


(A71 ) 


where  G(r)  =  J  pF(p)  sinh{^(p  -  r)}dp  . 


(A72) 


Now,  from  Eq.  (A62),  u(r)  =  Kje  Xr  +K  1okx  +  Pdr  +  AG(r),  and  using  the  transformation  P 
u/r,  the  general  solution  of  the  diffusion  Eq.  (A60)  is 


P(r,t)  =  Pd  +  ^^  +  K1— +  K2  — 
r  r  r 


(A73) 


The  integration  constants  Ki  and  K2  are  determined  by  applying  the  boundary 
conditions,  P(n,  t)  =  Pj(t)  and  P(r0,  t)  =  Pt(t)  to  Eq.  (A73).  We  obtain 


(e  )kj  +  (eXr‘  )k2  =-[(Pd  -Pi)ri  +70^)], 


(A74) 
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and  (e^r°  )Kj  +  (eIr°  )k2  =  -[(Pd  -  Pt)r0  +  AG(r0)] 

=  [(Pt  -Pd)r0-AG(r0)]. 

Solution  of  the  above  two  simultaneous  equations  in  Ki  and  K2  yields 

[(Pd  -Pi  X  +^G(r]  )]eXr°  +[(Pt  -pd  )r0  -AG(r0  )p 


(A75) 


hr: 


2  sinh(Ah) 


(A76) 


and 


k2  = 


[(Pd  -P;  )r,  +AG(r,  )]e-lr°  +  [(Pt  -Pd  )r0  -AG(r0 ) ]e~^ 
2  sinh(Ah) 


(A77) 


where  sinh(Ah)  =  - 


e+A,h  _  e“Xh 


,  and  h  =  r0  -  n. 


Substituting  the  above  expressions  for  Ki  and  K2  in  Eq.  (A73),  P(r,  t)  in  the  region 
n  <  r  <  r0,  is  given  by 


P(r,t)  =  Pd+^^ 


2  sinh(Ah)  r 
1  1 


1  1  Vd-P.  jr,  +  ^G(ri)jeA'(r°_r-)  +  {(Pt-  Pd)r0  - 

.{(Pd +  ^G(ri)je_A’(r°_r>  +  {(Pt-Pd)r0  - 


2  sinh(Ah)  r 


p  (  AG(r)  (Pd~  Pj)ri  +7G(rj)  sinh{A(r0  -r)} 
d  r  sinh(Ah)  r 

+  (Pt  -Pd)r0  ~AG(r0)  smh{A.(r  —  rt )} 
sinh(Ah)  r 


(A78) 


The  pressure  P(r,  t)  varies  with  time  due  to  changes  in  boundary  pressures  as  well  as 
any  changes  in  the  deviations  function  F(r)  with  time.  The  terms  involving  the  function 
G(r)  in  Eq.  (A78)  arise  from  the  non-uniform  component  F(r)  of  the  sink  pressure 
included  in  the  diffusion  Eq.  (A60).  These  terms  appear  in  combination  with  A  and 
therefore  vanish  if  A  =  0.  Thus, 


lim  P(r,t)  =  Pd  - 

X—>0 


(Pd-Pjfr 

h 


(rQ~r)  !  (PfPdK 

r  h 


(r  ~  ri  ) 
r 


Pd  -(Pd-Pi)ir 

h 


— 

r; 

-5--1 

+  (Pt-Pd)-r- 

1--L 

r 

h 

r 
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=  Pt  +(pd-pt)  -[(pt-pi)+(pd-pt)£ 


"*_r 

+  (Pt-Pd)v 

r 

h 

r 

T; 

r,  frA 

Pt  -(Pt -Pi )t- 

-5--1 

+  (Pd-Pt) 

l-t  -- 

h 

r 

h  [  r 

1*: 

r„ 

r„  -  r; 

i 

i 

=rT- 

-2--1 

r 

+  (Pd-Pt) 

1  0  1 
h 

o  1  1 


=  pt-(p,-pi)lL 


-1 


using  r0-  n  =  h, 


(A79) 


which  is  in  agreement  with  Eq.  (A10)  if  the  diffusion  region  is  not  perfused. 
Expressions  for  Pressure  Gradients  and  Gas  Fluxes 

We  see  from  Eq.  (A78)  that  the  gradient  dP/dr  involves  the  derivative  of  G(r),  which  we 
evaluate  using  Eq.  (A72): 


dG(r)  _  d 
dr  dr 

rr  d 


|  pF(p)  sinh{A,(p  -  r)}dp 


J  [PF(P)  sinh{^(p  -  r)}]dp  +  [pF(p)  sinhfMp  -  r)}]|p=r 


=  -  A.J  pF(p)  cosh{A,(p  -  r)}dp  +  0 
=  -  7  FT  ( r ) , 

where  H(r)  =  J  pF(p)cosh{^(p  -  r)dp  . 


(A80) 
(A81 ) 


Differentiating  P(r,  t)  given  by  Eq.  (A78)  with  respect  to  r,  and  using  Eq.  (A80),  we 
obtain 


<9P  _  ^G(r)  X  dG(r)  (Pd  -  Pj)^  +  ^G(ri) 

dr  r2  r  dr  sinh(^h) 


(Pt-Pd)r0-H3(r0) 


sinh(>uh) 


^cosh{^(r0  —  r)}  sinh{^(r0 -r)} 


?^cosh{^(r -r;)}  sinhf^r-r;)} 
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AG(r)  A2  H(r)  (Pd  -  fifr  +  AGO-) 

r2  r  sinh(Ah) 


Acosh{A(r0 -r)}  sinh{A(r0 -r)} 


,  (Pt-Pd)r0-XG(r0) 
sinh(Ah) 


Acosh{A(r  -r;)} 
r 


sinh{A(r-rj )} 


(A82) 


It  can  be  shown  that  this  gradient  expression  reduces  to  that  of  the  3RWT  models  given 
by  Eq.  (A1 1 )  if  A  =  0.  We  determine  the  gradients  gi  and  g0  by  evaluating  dP/dr  at  r  =  n 
and  r  =  r0,  respectively,  substituting  for  G(n),  G(r0),  H(n),  and  H(r0),  and  simplifying  the 
expressions  using  hyperbolic  function  identities  and  the  relationship  r0  =  (n  +  h): 


gi 


dP 

dr 


AG(r;)  A2H(ri) 


+ 


(Pd  -  Pi  )ri  +  AG(r;) 


sinh(Ah)  sinh(Ah) 


A  cosh(Ah)  sinh(Ah) 


"(Pt-PdK 

AG(r0)  " 

"a" 

sinh(Ah) 

sinh(Ah) 

_ri_ 

A2H(ri) 

ri 


+  (Pd 


cosh(Ah) 

sinh(Ah) 


1 

+  — 

ri. 


|  A2G(ri)cosh(Ah)  |  (Pt  -Pd)Ar0 
r;  sinh(Ah)  rs  sinh(Ah) 


A2G(r0) 

Tj  sinh(Ah) 


=  (Pd 


-Pi) 


cosh(Ah) 

sinh(Ah) 


1 

+  — 
ri . 


(Pt-Pd)Ar0  |  A2!, 

r;  sinh(Ah)  r;  sinh(Ah) 


(A83a) 


where  I;  =  -sinh(Ah)H(rj)  +  cosh(Ah)G(r;)  -  G(r0) ,  and 


8P 

AG(r0) 

A2H(rc)  , 

1 

1 

AG(rj)  ' 

"a" 

dr 

r2 

r=ro.  lO 

i 

ro 

sinh(Ah) 

sinh(Ah) 

_ro_ 

£ 

IP? 

1 

oT 

1 _ 1 

AG(r0)  ' 

A  cosh(Ah)  sinh(Ah) 

sinh(Ah) 

sinh(Ah) 

■.  d  \ 

A2H(r0)  |  (Pd  ~  Pj)Ar,  |  A2G(r,)  | 

r0  r0  sinh(Ah)  r0  sinh(Ah)  1 


cosh(Ah) 

sinh(Ah) 


A2G(r0)cosh(Ah) 
r0  sinh(Ah) 


(Pd-Pj)^i 

r0  sinh(Ah) 


+  (Pt~Pd) 


A 


cosh(Ah)  1 
sinh(Ah)  r0 


+  - 


A2Ir 


r0  sinh(Ah) 


where  I0  =  -sinh(Ah)H(r0)-cosh(Ah)G(r0)  +  G(ri) . 


(A83b) 
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The  integral  \\  is  evaluated  using  the  definitions  of  G(r)  and  H(r)  given  by  Eqs.  (A72)  and 
(A81 )  to  obtain: 


I;  =  -sinh(Fh)J  ‘  pF(p)cosh{F(p  -r;)}dp  +  cosh(Fh)J '  pF(p)sinh{F(p  -r;)}dp 
-J°pF(p)sinh{Mp-r0)}dp 


(A84) 


Replacing  r0  by  (n  +  h)  in  the  last  integral, 


F  =  - 


sinh(A,h)J '  pF(p)cosh{F(p  -r;)}dp  +  cosh(Fh)| '  pF(p)  sinh{F(p  -  r;)}dp 
-J  °  pF(p)  sinh {F( p  -r,  -  h)}dp 


(A85) 


Expanding  sinhfA.fp-r; -h)}  in  the  last  integral, 


I;  =  -  sinh(Fh)|  ‘  pF(p)cosh{F(p  -  r;)}dp  +  cosh(Fh)| '  pF(p)  sinh{F(p  -  r;)}dp 
-J  °  pF(p)[sinh{F(p  -  r^jcoshfTih)  -  cosh{F(p  -  ri)}sinh(Fh)]dp 


(A86) 


Rearranging  terms, 


I;  =sinh(Fh)  | '  pF(p)cosh{F(p-ri)}dp-| '  pF(p)cosh{F(p -r;)}dp 
-  cosh(Fh)|~ |  °  pF(p)  sinh{F(p  -  r;  )}dp  -  J" '  pF(p)  sinh{F(p  -  r;  )}dp 


=  sinh(Fh)Ic- cosh(Fh)Is , 


(A87) 


where  Ic  =  J°  pF(p)cosh{A.(p-ri)}dp  and  Is  =  J°pF(p)sinh{A,(p-ri)}dp  . 

Similarly,  the  definitions  of  G(r)  and  H(r)  given  by  Eqs.  (A72)  and  (A81 )  are  also  used  to 
evaluate  the  integral  l0: 

I0  =  -  sinh(A.h)|  °  pF(p)cosh{A.(p  -r0)}dp  -  cosh(A,h)|  °  pF(p)sinh{A,(p  -  r0)}dp 
+J‘pF(p)smh{A,(p-ri)}dp 


■  |  °  pF(p)[sinh(A.h)cosh{A.(p  -r0)}  +  cosh(A,h)sinh{A,(p  -  r0)}]dp  +  | '  pF(p)sinh{A,(p  -r;)}dp 
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J  °  pF(p)  sinh{7(p  -  r0  +  h)}dp  -  J '  pF(p)sinh{7(p  -  r,  )}dp 


(A88) 


Replacing  r0  by  (n  +  h)  in  the  first  integral, 


I 


O 


J  °  pF(p)sinh{Mp  -r;)}dp  - 1 '  pF(p)sinh{7(p  -  r;)}dp 


-L 


(A89) 


We  now  substitute  the  identities  of  li  and  l0  in  Eqs.  (A83a)  and  (A83b)  for  gi  and  g0  to 
obtain,  respectively: 


Si  =(Pd -Pi) 


X  coth(A,h)  + 


_1_ 

ri 


,  (Pt-Pd)^o 

Tj  sinh(7h) 


— — — — | [sinh(Xh)Ic  -  cosh(7h)ls] 
r;  sinn(A,h) 


=  (Pd 


+ 


'(Pt-Pd)fcr0 

<  - 

r;  sinh(Mi) 


+  — •  u/n,x[sinh^h)I^  “  cosh(A.h)Is  ]  +  (Pd  -  ^ )X coth(A,h) 
r;  smh(An) 


and 


§0 


(Pd  -Pj)foj 

r0  sinh(A,h) 


K  coth(A,h) 


^2IS 

r0  sinh(A,h) 


(A90a) 


(A90b) 


These  pressure  gradients  are  multiplied  by  the  corresponding  surface  areas  and  gas 
permeabilities  to  yield  the  gas  fluxes  fj  and  f0  at  the  inner  and  outer  boundaries  of  the 
diffusion  region: 


2 

D„ 

2 

fD  1 

fi  =«tDsgiAi  =atDsgi(4^rr)  = 

s 

_Db_ 

atDbgi(4™r)  = 

1 

z> 

Q 

_ 1 

and  f0  =  atDbg0A0  =  atDbg0(4ur02),  (A91b) 

where  fj0  is  the  gas  flux  at  inner  boundary  of  the  diffusion  region  coinciding  with  the 
bubble  surface.  This  flux  would  equal  fj  if  Ds  =  Db  with  no  differential  diffusivity,  and  it 
would  be  the  same  as  f0  in  the  absence  of  blood  flow  in  the  diffusion  region,  i.e. ,  in  the 
absence  of  sinks.  Recall  that  in  3RWT  models  with  an  unperfused  diffusion  layer  around 
the  gas  bubble  the  flux  in  the  diffusion  region  as  given  by  Eq.  (A17)  is  independent  of  r. 
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Equation  for  Rate  of  Change  of  Bubble  Radius 


Under  the  condition  of  zero  flux  at  the  outer  boundary  r0,  Eq.  (A90b)  is  solved  for  the 
integral  ls  with  gQ  =0  to  obtain 


X2 1 s  =  (Pd  -  Pj  )A,r;  +(Pd  -  Pt  )[sinh(Xh)  -  Xr0  cosh(Xh)] . 


(A92a) 


Because  the  flux  is  zero  at  the  outer  diffusion  boundary,  the  gradient  and  flux  at  the 
bubble  surface  must  be  independent  of  gas  tension  Pt  in  the  well-stirred  tissue  outside 
the  diffusion  region.  They  are  also  independent  of  the  thickness  h  in  this  large-scale 
model  where  the  diffusion  region  is  not  a  thin  layer  as  in  the  3RWT  model.  Nor  are  they 
dependent  on  gas  tension  deviations  due  to  F(r),  because  of  the  arbitrary  nature  of  F(r) 
allowed  to  accommodate  a  finite-sized  diffusion  region.  In  fact,  as  shown  later  [see  Eq. 
(A104)],  the  spatial  average  of  F(r)  decreases  with  increasing  diffusion  region  volume, 
reaching  zero  in  the  limit  at  infinite  volume  corresponding  to  the  2R  model.  Thus,  bubble 
growth  is  unaffected  by  tissue  gas  tension  Pt,  diffusion  region  thickness  h,  and  sink 
pressure  deviations  F(r).  As  all  dependence  of  the  bubble  surface  gradient  on  these 
factors  is  manifest  in  the  second  term  in  braces  on  the  right  side  of  Eq.  (A90a),  the 
gradient  is  rendered  independent  of  these  factors  by  setting  this  term  equal  to  zero.  An 
expression  for  the  gradient  consistent  with  that  for  the  2R  model  is  obtained  by  first 
rewriting  Eq.  (A90a)  in  equivalent  form; 


gi=(Pd-Pi) 


?i  +  — 
T: 


(Pt-Pd)K,  ,  'X 

+  \ - - - 1 - - - 

r;  sinh(A,h)  r;  sinh(A,h) 


and  setting  the  revised  term  in  braces  equal  to  zero: 


[sinh(A.h)Ic  -cosh(A.h)Is  +  ( Pd  —  P; )[A. coth(A-h)  —  A.] I ; 


(Pt-Pd)K 


-[sinh(A.h)Ic  -  cosh(7Ji)Is  +  (Pd  -  P, )[/. coth(Ah)  -  /.]  =  0 


rs  sinh(Xh)  r,  sinh(Xh) 

Substituting  A.2IS  from  Eq.  (A92a)  into  Eq.  (A92b)  and  solving  for  ^2IC  yields 

/.2IC  =  (Pd  -  P,  )/.r,  +  (Pd  -  Pt  )[cosh(Xh)  -  /.r0  sinh(A-h)] . 


(A92b) 


(A92c) 


The  complementary  Eqs.  (A92a)  and  (A92c)  characterize  F(r),  and  reduce  Eqs.  (A90a) 
and  (A90b)  to 


0P 

dr 


1 

\  11 

II 

II 

D. 

1 

SP 

1 

+ 

1 _ 

,  and 


(A93a) 
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r=rf 


o ,  respectively. 


(A93b) 


8P_ 

dr 


Eq.  (A93a)  is  substituted  for  gi  in  Eq.  (A4)  to  obtain  the  expression  for  the  rate  of 
change  of  bubble  radius: 


dr^ 

dt 


a'tDs(Pd-Pi) 


X  +  - 


*~i  ^amb 

3  dt 


amb 


-P: 


idg 


4y 

+  — + 

3rj 


8n 

3 


Mrf 


(A94) 


Note  that  the  bubble  surface  diffusivity  Ds  may  be  different  from  the  bulk  diffusivity  Db 
due  to  altered  diffusion  characteristics  of  the  bubble  surface. 


Average  Tissue  Gas  Tension 

The  equation  for  determining  the  gas  tension  in  the  diffusion  region  is  based  on  the 
following  expression  for  mass  balance  in  this  region: 

r  r 

at  ~|p(pJt)47tp2dp=ab |[(Qt  +Qn(p,t)][Pa  -P(p> t)]4np2dp-fi ,  (A95) 

ri  ri 


where  the  diffusion  region  gas  content  on  the  left  of  the  equation  is  from  Eq.  (A32),  Qtis 
the  uniform  blood  flow  per  unit  tissue  volume  equal  to  the  average  blood  flow  per  unit 
tissue  volume,  Qn(p,t)  is  the  position-dependent  non-uniform  component  of  blood  flow 
per  unit  tissue  volume  representing  perfusion  heterogeneity,  which  can  vary  with  time, 
and  fi  is  the  flux  into  or  out  of  the  bubble.  Note  that  Eq.  (A95)  holds  regardless  of  the 
shape  of  the  diffusion  region. 

Expressing  the  flux  fi  as  the  rate  of  change  of  bubble  gas  content  from  the  left  side  of 
Eq.  (A1),  Eq.  (A95)  becomes 

r  r 

“i^jpCMMTtrdp^hQi  j[pa  -P(p,t)]4Itp2dp-^rA(pivi)  +  e(t)  (A96) 

r  r 

Ai  Ai 


where  e(t)  =  ab  jQn(p,t)[Pa  -P(p,t)]47ip2dp  is  the  net  gas  content  carried  by  blood  into  or 

ri 


out  of  the  diffusion  region  by  the  heterogeneous  perfusion  components. 
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The  integral  on  the  left  side  of  Eq.  (A96)  evaluates  to  PdVd,  and  that  on  the  right  side  to 
(pa  -Pd)vd,  where  Vd  is  the  volume  of  the  diffusion  region.  Assuming  s(t)  =  0, 
substituting  the  evaluated  expressions  for  the  integrals  in  Eq.  (A96),  and  dividing  both 
sides  by  at,  yields 


Apjyd)= 

at 


(Pa"Pd)Vd 

T 


J_d_ 

a't  dt 


(Pi  v;). 


(A97) 


Equation  (A97)  involves  the  diffusion  region  volume  Vd,  which  may  change  during 
bubble  growth.  However,  we  choose  to  treat  Vd  as  a  constant  and  attribute  any  effects 
of  its  changes  to  changes  in  F(r).  Eq.  (A97)  for  the  average  diffusion  region  gas  tension 
then  simplifies  to 


dPd  _Pa~Pd 

dt  x 


1  d 
a'tVd  dt 


(Pi  Vi). 


(A98) 


For  large  diffusion  region  volumes,  the  second  term  on  the  right  side  of  Eq.  (A98)  drops 
out,  the  average  gas  tension  Pd  is  the  same  as  the  uniform  gas  tension  in  the  well- 
stirred  tissue  far  away  from  the  bubble,  and  the  3RUT  model  reduces  to  the  2R  model. 
Note,  however,  that  until  this  limit  is  reached,  Pd  will  generally  not  equal  the  gas  tension 
Pt  in  the  well-stirred  region,  which  varies  independently  of  gas  fluxes  into  or  out  of  the 
bubble  according  to  Eq.  (A54). 


Minimum  Diffusion  Region  Volume 

The  choice  of  constant  Vd  must  be  made  within  certain  constraints,  constraints  that 
follow  from  the  expression  for  the  gas  content  of  the  diffusion  region  and  the 
relationship  between  Vd  and  the  sink  pressure  deviations  function  F(r).  Substituting  the 
identity  for  the  blood-tissue  gas  exchange  time  constant,  the  diffusion  region  gas 
content  Ud  given  by  Eq.  (A32)  becomes 

Ud  =  at  J  °  P(p,t)4np2dp  =  xabQ,  J  °  P(p,t)4rtp2dp  .  (A99) 

The  integral  on  the  right  side  of  Eq.  (A99)  can  be  evaluated  directly  using  Eq.  (A78)  for 
P(r,  t),  but  the  derivation  is  tedious.  We  can  proceed  more  simply  by  noting  that  the  flux 
from  diffusion  region  sinks  must  equal  the  difference  (f0  -  f0)  between  the  gas  fluxes  at 
the  two  boundaries  of  the  region,  where  f0  is  the  flux  at  the  outer  boundary  and  fj0  is  the 
flux  at  the  inner  boundary  coinciding  with  the  bubble  surface.  Recall  that  fi0  may  differ 
from  the  flux  f  into  the  bubble  due  to  differential  diffusivity  with  Ds  *  Db  [see  Eq.  (A36)]. 
We  therefore  have 
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(A100) 


abQt  f0[p(p5t)-Ps(p)]4np2dp  -  f0  -fi0  =  f0 


Dt 

Dc 


-f; 


where  we  have  related  fi0  to  fi  using  Eq.  (A91a).  With  f0  =  0,  Eq.  (A100)  indicates  that 
when  the  bubble  is  expanding,  fi  >  0,  and  hence  the  average  sink  pressure  Ps(r)  is 

greater  than  the  spatial  average  of  P(r,t) ,  which  equals  Pd  by  definition,  to  serve  as  a 
source  of  gas  for  the  bubble;  the  opposite  is  the  case  when  the  bubble  is  shrinking. 

Now,  from  Eqs.  (A99)  and  (A100), 


Ud  =  Tab 


QtJr°ps(p)47TP2dP  +x 


fo- 


Dt 

Dc 


f; 


=  at|r°  Ps  (p)47Tp2dp  -T 


D| 

D, 


-f.  -f 

1  xO 


.  (A101 ) 


Substituting  for  Ps(r)  from  Eq.  (A58)  and  evaluating  the  integral, 


ud=  atVd(Pd  +  F)-x 


Dt 

Ds 


f.  -f 


(A102) 


where  F=—  fr°  F(p)47ip2dp  is  the  spatial  average  of  F(r),  and  (pd  +f)  is  the  average  sink 

Vd  *r 

pressure  Ps(r). 


The  foregoing  derivation  of  Eq.  (A102)  does  not  depend  on  the  shape  of  the  diffusion 
region  nor  on  a  spherically  symmetric  diffusion  process  with  radial  gas  fluxes.  We  use 
this  fact  in  a  later  section  to  extend  the  3RUT  model  to  conditions  that  involve  multiple 
bubbles  (3RUT-MB  model). 

Since  the  average  gas  tension  Pd  in  the  diffusion  region  is  given  by  ud  =  atPdVd ,  we 
obtain  the  following  from  Eq.  (A102): 


atVd 


Dt 

D„ 


-f:  -f„ 


(A103) 


Letting  f0  =  0,  and  expressing  f  as  the  rate  of  change  of  bubble  gas  content  from  the  left 
side  of  Eq.  (A1 ),  Eq.  (A1 03)  becomes 


F  = 


Dh  d 


a,(Vd  Ds  dt 


(piVi). 


(A104) 


Equation  (A104)  shows  that  the  spatial  average  gas  tension  deviation  F  remains 
positive  during  bubble  expansion  due  to  increasing  bubble  gas  content,  and  turns 
negative  during  bubble  contraction,  eventually  reaching  zero  when  the  bubble  has 
completely  resolved.  It  defines  a  hypothetical  ‘point’  source  or  sink  that  accounts  for  gas 
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diffusion  from  all  distributed  sources  or  sinks  within  the  diffusion  region.  Note  that  f 
gets  progressively  smaller  with  larger  values  of  Vd,  but  the  product  |F|vd  is  finite  at  any 
time  during  the  lifetime  of  the  bubble,  regardless  of  the  size  of  Vd. 

We  obtain  the  following  expression  for  F  by  using  Eqs.  (A1 ),  (A93a),  and  (A2)  in  that 
order  on  the  right  side  of  Eq.  (A104),  and  expressing  A,  in  terms  of  n: 


—  47r:ri2TDb 

",  1  " 

p  _  1  0 

X  +  — 

Pd  - 

vd 

ri_ 

^  amb 


2y  47iMri?  X 

pidg  +  — +  — 

'  j 


(A105) 


As  the  average  deviation  f  varies  to  satisfy  Eq.  (A105)  with  constant  Vd,  Vd  must 
remain  larger  than  a  certain  minimum  value  to  keep  the  average  sink  pressure 
ps  =  pd  +  f  realistic.  Values  of  Vd  below  this  minimum  will  result  in  values  of  |f|  large 

enough  to  render  ps  negative  during  bubble  resolution  when  F<0.  Also,  depending  on 
Vd,  Psmay  become  unrealistically  large  during  bubble  growth  when  F>0.  Thus,  Eq. 
(A105)  indicates  that  for  given  values  of  tissue  time  constant,  diffusivity,  and  surface 
tension,  there  is  a  minimum  diffusion  region  volume  required  for  bubble  evolution. 
Equivalently,  for  a  given  volume  Vd,  bubble  evolution  can  occur  only  if  the  diffusivity  is 
below  a  certain  maximum  value,  which  clearly  implies  diffusion  limitation. 

The  lower  bound  on  the  diffusion  region  volume,  which  we  will  denote  by  Vd(min),  can  be 
determined  for  any  given  decompression  profile  by  varying  Vd  with  all  other  model 
parameters  held  constant  and  imposing  the  constraint  that  (pd  +f)  remain  non-negative 
at  all  instants  during  the  lifetime  of  the  bubble.  A  given  minimum  value  for  Vd  chosen  at 
one  point  in  a  decompression  profile  may  require  adjustment  at  subsequent  points  in 
the  profile  where  changes  in  breathing  gas  composition,  pressure,  etc.,  occur.  Even  so, 
Vd  is  piece-wise  constant,  so  that  Eq.  (A98)  will  still  apply  after  each  adjustment. 

Equations  (A2),  (A94),  and  (A98)  together  with  \l\  -  (47tn3/3)  constitute  the  3RUT  model 
for  describing  the  evolution  of  a  single  bubble  in  a  tissue  with  finite  volume.  The 
adoption  of  perfusion  heterogeneity  in  the  diffusion  region  overcomes  the  requirement 
for  an  infinite  tissue  volume  in  the  2R  model,  but  the  assumptions  of  a  spherically 
symmetric  diffusion  process  and  constant  Vd  preclude  any  interactions  between 
different  bubbles  in  a  given  tissue  volume:  Each  bubble  evolves  within  its  own 
independent  diffusion  region  with  zero  flux  at  the  outer  boundary  that  everywhere 
adjoins  the  well-stirred  region.  The  model  is  therefore  only  a  single  bubble  model. 

The  3RUT  model  requires  simultaneous  solution  of  only  two  ordinary  differential 
equations.  The  model  parameters  to  be  specified  are:  y,  t,  a;,  Vd,  Ds,  Db  and  the  initial 

critical  radius  rc  (bubble  nucleus  size).  The  sink  coefficient,  X,  is  given  by^/l/xDb  ,  using 
the  definitions  of  X  and  t  under  Eqs.  (A46)  and  (A16),  respectively.  The  parameter  a't 
always  appears  in  combination  with  either  Ds  [Eq.  (A94)]  or  Vd  [Eq.  (A98)]  and  therefore, 
need  not  be  specified  independently.  Also,  we  will  have  one  less  free  parameter  if  we 
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assume  Ds  =  Db.  Numerical  implementation  of  Eq.  (A98)  is  the  same  as  that  for  Eq. 
(A42)  in  the  3RWT-VTDD,  but  without  the  necessity  for  any  approximation  since  the 
diffusion  region  volume  does  not  vary  continuously  with  bubble  evolution.  Note  that  Eq. 
(A105)  for  the  average  deviation  f  is  not  involved  in  determining  the  bubble  radius,  and 
is  needed  only  to  ensure  that  Vd  remains  sufficiently  large  to  keep  Ps  within  realistic 
limits  during  the  lifetime  of  the  bubble. 


V.  Three-Region  Unstirred  Tissue,  Multiple  Bubble  Model  (3RUT-MB) 

The  assumption  of  spherically  symmetric  diffusion  around  the  bubble  is  relaxed  in  the 
3RUT-MB  model,  as  is  the  assumption  of  constant  diffusion  region  volume  Vd.  We  begin 
by  considering  the  general  form  of  the  diffusion  equation  in  spherical  coordinates  with 
an  added  sink  term,  but  without  the  explicit  time-dependent  term,  as  allowed  by  the 
quasi-static  approximation: 

V2P(r,e^)=X2[P(r,0,^)-Ps] ,  (A106) 

where  P  is  tissue  gas  tension  at  any  radial  distance  r  from  the  center  of  the  bubble,  0 
and  \| /  are  angular  (polar  and  azimuth)  variables  of  the  spherical  coordinate  system,  A,  is 
a  constant  as  defined  below  Eq.  (A46)  for  all  unstirred  tissue  models,  and  Ps  is  the  sink 
pressure  given  by 


Ps(r,0,f)  =  Pd+F(r,e,f>.  (A107) 

As  in  Eq.  (A58)  for  the  3RUT  model,  Pd  is  the  location-independent  uniform  component 

of  Ps,  but  F(r,  0,  <\>)  is  now  any  function  of  r,  0,  and  (j).  The  uniform  component  is  the 
spatial  average  gas  tension  in  the  diffusion  region,6  which  may  be  taken  to  be  the  same 
as  the  average  gas  tension  of  the  entire  tissue  with  multiple  bubbles  (as  indicated  later). 
The  function  F(r,  0,  ()))  represents  the  non-uniform  component  of  Ps  arising  from 
deviations  in  sink  pressure  from  its  uniform  component  due  to  non-uniform  blood  flow. 
As  in  the  3RUT  model,  F(r,  0,  c)))  and  Ps(r,  0,  ()>)  are  functions  of  time,  but  we  omit  explicit 
indication  of  the  time-dependencies  to  keep  the  notation  simple. 

It  is  not  possible  to  solve  Eq.  (A106)  analytically  for  gas  tension  P  with  any  arbitrary 
F(r,  0,  §).  However,  bubble  growth  and  resolution  is  governed  only  by  the  gas  flux  into 
or  out  of  the  bubble,  and  hence  by  the  radial  components  of  the  gas  tension  gradient  at 
the  bubble  surface,  regardless  of  the  processes  beyond  the  bubble  surface  that  give 
rise  to  those  components.  We  can  therefore  exploit  the  arbitrary  nature  of  F(r,  0,  to 
specify  its  form  as  required  to  apply  the  solution  for  the  gradient  at  the  bubble  surface 
given  by  the  radially  symmetric  solution  of  the  diffusion  equation.  As  F(r,  0,  §)  varies 
with  direction  (0,  ()>),  the  outer  boundary  of  the  diffusion  region  will  also  vary  and  deviate 


cThe  "bar"  notation  to  indicate  spatial  average  was  omitted  in  the  original  description  of 
this  model  in  reference  6. 
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from  radial  symmetry.  The  applicable  solution  of  the  diffusion  equation  is  consequently 
that  for  P  in  a  concentric  near  diffusion  region  of  arbitrarily  small  thickness  immediately 
surrounding  the  bubble  where  radial  symmetry  can  be  assumed.  With  elimination  of  the 
requirement  for  radial  symmetry  in  the  remainder  of  the  diffusion  region,  the  3RUT 
model  limitation  that  the  outer  boundary  of  the  diffusion  region  must  everywhere  adjoin 
the  well-stirred  outer-most  tissue  region  is  overcome. 

With  inner  boundary  of  the  near  region  at  the  bubble  surface  n  and  outer  boundary  at 
ron,  the  integration  constants  for  the  radially  symmetric  solution  of  the  diffusion  equation 
in  Eq.  (A73)  are  evaluated  at  boundary  conditions  P(n,  t)  =  Pj(t)  and  P(ron,  t)  =  Px(t)  to 
obtain  the  following  expression  for  P(r,  t)  valid  for  n  <  r  <  ron: 

p(r  t)  =  p  +  ^G(r)  _  (Pd-Pj)ri  +7G(rj)  sinh{Mron  - r)} 

1  ^  d  r  sinh(/U|)  r 

(A108) 

+  (Px  -PdKn  -^G(v)  sinh{A.(r  —  rt )} 
sinh(^h)  r 


where  Pi  is  the  bubble  gas  pressure,  h  =  (ron  -  n)  is  the  thickness  of  the  near  diffusion 
region,  Px  is  the  gas  tension  at  r  =  ron,  and  G(r)  is  as  defined  in  Eq.  (A72).  Differentiating 
Eq.  (A1 08)  with  respect  to  r  yields 


dP 

dr 


XG(r) 


72H(r)  |  (Pd  -PQrj  +kG(ri) 
r  sinh(/Ui) 


sinh{7(ron 


r)}  +  hr  cosh  {7(ron  -r)} 


(A109) 


(pd  -pxKn  +^G(ron) 


sinh(7h) 


sinhfA.fr-r;)  -  Xr  cosh{A.(r  -r;)} 


where  H(r)  =  ~  — dG^  ,  as  defined  by  Eq.  (A80).  Eq.  (A109)  is  evaluated  at  r=n  and  r=ron 

X  dr 

to  obtain  the  following  gradients  at  the  inner  and  outer  boundaries,  respectively,  of  the 
near  region: 


ap 

dr 


r=ri 


=  (Pd-Pi) 


X  coth(7h)  + 


ri 


(pd-px)7ron 
Tj  sinh(/,h) 


—  [ls  cosh(7h)  -  Ic  sinh(7h)] 

r;  smh(A,h) 


(AllOa) 


and 
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Ic  =  J  PFn  (P)  cosh{A,(p  -  r,  )}dp 

rj 

with  Fn(p)  representing  the  deviations  function  F(r,  0,  (j>)  in  the  near  region.  The 
gradients  given  by  Eqs.  (A1 10a)  and  (A1 10b)  can  be  satisfied  by  appropriate 
specification  of  the  deviations  function  Fn(r),  or  equivalently,  by  specifying  A.2IC  and  A,2IS 
at  every  point  in  the  diffusion  region.  Denoting  the  gradient  at  ron  by  gx,  Eq.  (A1 10b)  is 
solved  for  the  integral  ls  to  obtain: 

^2IS  =  (Pd  -Pj )Xr;+(Pd  -  Px  )[sinh(>Ji)  -  Xron  cosh(/.h)]-gxron  sinh(Ah) .  (A1 11a) 


Spherically  Symmetric  Fn  in  the  Near  Region 

If  the  function  Fn(r,  0,  (j))  is  assumed  to  be  spherically  symmetric,  the  gradient  and  flux  at 
the  bubble  surface  are  independent  of  0  and  (j).  The  gradient  and  flux  at  the  bubble 
surface  are  also  independent  of  h  and  the  gas  tension  Px  at  ron,  both  of  which  are 
arbitrary.  Nor  are  they  dependent  on  gas  tension  deviations  due  to  Fn(r),  because  of  the 
arbitrary  nature  of  the  deviations  function  allowed  to  accommodate  a  finite-sized 
diffusion  region.  Thus,  bubble  growth  or  resolution  is  unaffected  by  gas  tension  Px, 
diffusion  region  thickness  h,  and  sink  pressure  deviations  Fn(r).  This  is  ensured  by 
specifying  the  second  term  in  braces  on  the  right  side  of  Eq.  (A1 10a)  as 

-  -  •  Trrr  ~ TTTTT^  cosh(^)-Ic  sinhfXh)]  =  (Pd  -  P, )[/.  coth(Xh)  -  /,] .  (A1 11b) 

r;  sinh(Ah)  r;  smh(Ah) 

As  in  the  derivation  of  Eq.  (A92b)  for  the  3RUT  model,  the  factor  {Xcoth(Xh)-X)  on  the 
right  side  of  Eq.  (A1 1 1b)  eliminates  h  dependence  in  the  expression  for  the  bubble 
surface  gradient,  here  given  by  Eq.  (A1 10a),  and  leads  to  an  expression  for  that 
gradient  consistent  with  the  corresponding  expression  in  the  2R  model.  Substituting  A,2IS 
from  Eq.  (A1 1 1  a)  into  Eq.  (A1 1 1  b)  and  solving  for  X2\c,  we  obtain 

^Ie=(Pd  -Pi)^i+(Pd-px)[cosh(^h)  -^ron  sinh(Mi)]-gxron  cosh(Mi) .  (A1 1 1  c) 
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The  complementary  Eqs.  (A1 1 1  a)  and  (A1 1 1  c)  characterize  Fn(r)  in  the  near  region,  and 
reduce  Eqs.  (A1 10a)  and  (A1 10b),  respectively,  to 


SP 

dr 


r  ii 

=(pd-pi) 

r=rj 

i 

+ 

l _ 

and 


dP 

dr 


§X  ' 


(All  2a) 


(All  2b) 


Eq.  (A1 12a)  is  equivalent  to  the  3RUT  model  solution  for  the  gradient  at  the  bubble 
surface.  The  inclusion  of  the  term  containing  gx  in  Eqs.  (A1 1 1  a)  and  (A1 1 1  c)  does  not 
alter  the  model  from  the  3RUT  model  for  a  single  bubble:  The  3RUT  condition  of  gx  =  0 
is  simply  a  particular  case  of  the  current  more  general  formulation. 


Spherically  Asymmetric  Fn  in  the  Near  Region 


We  now  relax  the  assumption  of  spherically  symmetric  Fn.  Substituting  A.2IS  from  Eq. 
(A1 11a)  and  simplifying  using  the  hyperbolic  identity  {cosh(A.h)}2  - (sinh(A.h)}2  =1 ,  Eq. 
(A1 10a)  becomes 


dP 

a7 


r=ij 


=  (Pd-Pi) 


(Px  -  Pd)lXron  sinh(Tih)  -  cosh(/,h )]  |  gxron  cosh(lh) 

ri  ri  J 


Eq.  (A1 13)  expresses  the  radial  component  of  the  gradient  at  any  point  on  the  bubble 
surface,  regardless  of  radial  direction  (9,  (j)),  implying  that  it  is  independent  of  the 
quantities  h,  Px,  gx,  and  lc  associated  with  the  point  along  any  given  radius  at  the  outer 
boundary  ron.  Thus,  the  term  in  braces  on  the  right  side  of  Eq.  (A1 13)  must  be 
dependent  only  on  0  and  (j).  We  therefore  define 


r(e,  <t») = 


(Px  -pd)[^ron  sinh(Xh)-cosh(Xh)]  |  gxron  cosh(A.h) 


(A114) 


where  r  is  a  function  only  of  0  and  (j),  so  that  Eq.  (A1 13)  becomes 


dP 

dr 


r=ri 


=  (Pd  “Pi) 


+  r(e,<|>). 


(All  5) 


‘  Note  that  substitution  of  A.2IC  from  Eq.  (A1 11c)  appropriately  reduces  Eq.  (A1 13)  to 
Eq.  (All 2a). 
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All  influences  of  diffusion  and  perfusion  on  the  gradient  beyond  n  in  the  (0,  4>)  direction 
are  now  represented  by  a  single  function  r(0,  c|>) . 

The  total  flux  at  the  bubble  surface  is  given  by 

fi=JKs§ids>  (All  6) 

S; 


where  Si  denotes  the  bubble  surface,  ds  is  the  elemental  surface  area,  and  ks  is  the 
gas  permeability  of  the  bubble  surface.  Substituting  the  gradient  from  Eq.  (A1 1 5)  and 
completing  the  integration,  we  obtain 


fi=Ks(Pd-Pi) 


4nr. 


•2  +  jKsr(e,(j))ds, 


(All  7) 


where  ics  is  the  average  gas  permeability  of  the  bubble  surface.  Without  loss  of 
generality,  we  express  the  integral  of  Ksr(e,<|>)  in  terms  of  quantities  in  the  first  term  on 
the  right  side  of  Eq.  (A1 1 7)  together  with  an  unknown  parameter  A: 

|  Ksr(e,  4>)ds  =  Aks  (Pd  -  Pj  )47ir2  .  (A1 1 8) 

Si 


Substituting  Eq.  (A1 1 8)  into  Eq.  (A1 1 7)  yields 


fj  =  atDs(Pd  -P;) 


A  +  - 
r; 


4 nr2  =  atDsg;Ai , 


(All  9) 


where  g;  =(Pd-P;) 


is  the  gradient  at  the  bubble  surface  and  ks  has  been 


replaced  by  the  product  of  the  average  gas  solubility  in  tissue  and  the  average  bubble 
surface  diffusivity;  ks  =  ajDs  =  atDs ;  for  convenience. 


The  gradient  given  by  Eq.  (A1 19)  is  the  average  over  the  bubble  surface  with  deviations 
from  this  average  due  to  differences  in  Fn(r,  0,  (j>)  canceling  each  other  out.  The  only 
assumption  remaining  besides  the  quasi-static  approximation  is  that  the  gas  molecules 
enter  and  leave  the  bubble  surface  in  the  normal  direction  perpendicular  to  the  surface, 
leading  to  an  average  bubble  surface  gradient  consistent  with  those  in  the  2R  and 
3RUT  models  with  spherically  symmetric  diffusion;  the  former  with  uniform  blood  flow  in 
large  tissue  volumes  and  the  latter  with  heterogeneous  blood  flow  in  finite  tissue 
volumes.  The  parameter  A  in  the  present  representation,  however,  can  validly  assume  a 
far  wider  range  of  values  than  the  corresponding  X  parameter  ( =  ^/abQt  /atDb  =  ^/l / xD b  ) 
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in  these  other  two  models.  Because  A  includes  influences  of  contributions  to  the  bubble 
surface  gradient  arising  from  the  radial  components  of  gas  fluxes  with  components 
tangential  to  the  bubble  surface,  it  may  assume  larger  values  than  the  applicable  X.  On 
the  other  hand,  A  may  assume  smaller  values  than  the  applicable  X  in  multiple  bubble 
systems  due  to  the  effects  of  bubble-bubble  interactions.  With  constant  values  for  the 
time  constant  x  and  diffusivity  Db,  X  is  a  constant,  but  not  A,  which  is  generally  time 
dependent,  reflecting  changes  in  bubble-bubble  interactions  caused  by  changes  in 
bubble  size  and  number  density. 


Equation  for  Rate  of  Change  of  Bubble  Radius 

If  the  diffusion  process  around  the  bubble  is  assumed  to  be  spherically  symmetric, 
substitution  of  the  gradient  from  Eq.  (112a)  for  gi  into  Eq.  (A4)  yields  Eq.  (A94),  the 
3RUT  model  expression  for  the  rate  of  change  of  bubble  radius. 


For  the  more  general  case  in  which  the  assumption  of  spherical  symmetry  is  relaxed,  A 
is  assumed  constant  so  that  the  expression  for  the  rate  of  change  of  bubble  radius  is 
consistent  with  the  corresponding  expressions  in  the  other  models.  Thus,  the  equation 
for  the  rate  of  change  of  bubble  radius  is  obtained  by  substituting  the  gradient  given  by 
Eq.  (A1 19)  into  Eq.  (A4)  to  yield 


dr, 

dt 
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(A120) 


Eq.  (A1 20)  is  identical  to  Eq.  (A94)  with  substitution  of  A  for  X.  As  noted  earlier,  Eq. 
(A94)  is  applicable  to  only  a  single  bubble  with  symmetric  diffusion  around  it  and  X 
determined  by  the  time  constant  x  and  diffusivity  Db.  Eq.  (A120),  on  the  other  hand,  is 
applicable  to  multiple  bubbles  in  the  tissue  and  must  be  solved  for  each  bubble  size  if 
the  bubbles  are  of  different  sizes.  We  will  denote  the  bubble  radii  by  (n)n,  n  =  1, 2,  ...Ns, 
where  Ns  is  the  number  of  different  bubble  sizes. 

Note  that  if  A  is  constant,  it  may  be  used  as  a  scaling  parameter  to  scale  bubble  radius 
and  other  model  parameters,  yielding  multiple  realizations  of  n(t)  with  a  single  solution 
of  Eq.  (A120).6’14 


Equation  for  Tissue  Gas  Tension 

While  we  require  the  near  region  to  be  spherically  symmetric  in  order  to  solve  the 
diffusion  equation  in  a  simple  fashion,  the  larger  diffusion  region  may  be  of  any  shape 
because  F  is  not  generally  spherically  symmetric.  The  total  thickness  of  the  diffusion 
region  may  therefore  vary  with  direction  as  well  as  with  time.  In  any  given  radial 
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direction,  diffusive  exchange  of  gas  between  bubble  and  tissue  proceeds  until  zero 
gradient  is  attained,  which  will  occur  where  the  diffusion  region  ends  in  contact  with  the 
well-stirred  outer  region  or  where  the  diffusion  region  abuts  that  of  an  adjoining  bubble 
at  a  non-spherical  interface.  At  such  regions  of  abutment,  the  direction  of  diffusive  gas 
exchange  changes  towards  the  other  bubble  allowing  the  bubbles  to  interact  with  one 
another  by  depriving  one  another  of  diffusion  region  space.  Accommodation  of  such 
interactions,  following  from  the  asymmetry  of  the  deviations  function,  is  the  feature  of 
this  model  that  distinguishes  it  from  the  3RUT  model  with  its  radially  symmetric 
deviations  function.  However,  the  exact  shape  of  the  diffusion  region  is  of  no  concern 
for  determining  mass  balance.  With  zero  flux  at  the  outer  boundary,  the  following  3RUT- 
MB  model  analog  of  Eq.  (A95)  applies  for  mass  balance  in  the  diffusion  region  around 
each  bubble: 

fp(z»t)dv  =  ab  j[(Qt  +  Qn(z>t)IPa  -P(z,t)]dv-fj,  (A121 ) 

dt  vd  vd 


where  Vd  is  the  volume  of  the  diffusion  region,  z  represents  the  coordinates  (r,  0,  and  ijj) 
of  any  point  within  the  diffusion  region,  dv  is  the  elemental  volume  of  integration,  Qt  is 

the  average  blood  flow  per  unit  tissue  volume,  and  Qn(z,t)  is  the  position-dependent 
non-uniform  component  of  blood  flow  per  unit  tissue  volume  representing  perfusion 
heterogeneity,  which  can  vary  with  time.  As  in  the  3RUT  model,  we  assume 
equilibration  of  tissue  and  venous  gas  tensions. 

Expressing  the  flux  as  the  rate  of  change  of  bubble  gas  content  from  the  left  side  of 
Eq.  (A1 )  and  using  the  suffix  k  to  denote  the  kth  bubble,  Eq.  (A1 21 )  becomes 


at 


Jpk (z,t)dv  =  abQt  J[Pa 

(vd)k  (vd)k 


-Pk(z,t)]dv+  JabQnk(z,t)[Pa 
(vd)k 


Pk  (z,  t)]dv 


J_A 

RT  dt 


(P,Vi)k. 


(A122) 


We  sum  both  sides  of  Eq.  (A122)  over  k  to  include  all  bubble  regions  (k  =  1  through  N 
where  N  is  the  total  number  of  bubbles  with  Ns  different  sizes)  and  ‘bubble-free’  regions 
(k  =  0),  for  which  the  last  term  on  the  right  side  is  zero,  to  obtain  the  following 
expression  for  mass  balance  in  the  entire  tissue: 


j  N  N 

at“X  JPk(z,t)dv  =  abQt^  J[Pa 
k=0(vd)k  k=°(vd)k 


1  w  H 

-Pk(z,t)]dv-  — ^  — (PiVi)k  +e(t) , 
RT  dt 


(A123) 


N 

where  e(t)  =  ^  JabQn  k(z,t)[pa -Pk(z,t)]dv  (A124) 

k=0(vd)k 
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is  the  net  gas  content  carried  by  blood  into  or  out  of  the  tissue  by  the  heterogeneous 
perfusion  components. 


With  pt  defined  as  the  spatial  average  gas  tension  in  the  entire  tissue,  including  the  net 
contribution  s(t)  from  the  heterogeneous  perfusion  components,  the  sum  of  integrals  on 
the  left  side  of  Eq.  (A123)  evaluates  to  PtVt,  and  the  sum  of  integrals  on  the  right  side  to 
(pa  -Pt)vt,  where  Vt  is  the  total  tissue  volume.  The  latter  is  fixed  and  equal  to  the  sum  of 
the  volumes  of  the  diffusion  regions  around  all  bubbles  and  the  volumes  of  the  well- 
stirred  bubble-free  regions.  Substituting  the  evaluated  expressions  for  the  integral  sums 
in  Eq.  (A123)  and  dividing  both  sides  by  atVt,  yields 


dPt 

dt 


(Pa-Pt) 


N 


atVt 


k=l 


s(t) 

atVt 


(A125) 


Note  that  (Vd)k  does  not  appear  in  Eq.  (A125)  as  the  diffusion  region  volume  of  each 
bubble  can  change  indeterminately  during  the  lifetime  of  the  bubble.  In  the  course  of 
such  changes,  the  diffusion  region  exchanges  volume  with  the  well-stirred  region  or  with 
the  diffusion  region  of  one  or  more  adjoining  bubbles.  Consideration  of  the  events  that 
must  cause  such  exchanges  compels  identification  of  the  average  diffusion  region  gas 
tension  Pd  as  the  gas  tension  Pt  in  the  well-stirred  region,  from  which  it  follows  that  Pd 
must  be  the  same  for  all  bubbles  at  any  given  time,  change  uniformly  across  the  entire 
tissue,  and  always  equal  the  prevailing  overall  average  tissue  gas  tension  pt . 

For  conceptual  consistency,  changes  of  bubble  diffusion  region  volume  cannot  involve 
any  discontinuity  in  sink  pressure  or  blood  flow  that  might  imply  existence  of  distinct 
anatomical  structure  at  the  outer  boundary  of  the  region,  or  any  global  changes  in 
properties  of  the  region,  as  the  volume  change  may  be  driven  by  localized  events  at  the 
periphery  of  the  region. 

If  a  gas  flux  to  or  from  a  bubble  develops  in  an  adjoining  volume  of  well-stirred  tissue, 
the  diffusion  region  of  the  bubble  expands  by  acquisition  of  volume  from  the  well-stirred 
region.  If  Pt  does  not  equal  Pd ,  such  a  flux  can  develop  only  with  a  discontinuous 
change  in  the  non-uniform  component  of  the  sink  pressure  F  in  the  well-stirred  tissue 
adjoining  the  original  boundary  between  the  diffusion  and  well-stirred  regions. 

Moreover,  the  tension  Pd  in  the  diffusion  region,  and  therefore  the  sink  pressure 
Ps  =pd  +f  ,  must  change  throughout  the  entire  diffusion  region.  Both  of  these 
requirements  are  unrealistic.  On  the  other  hand,  if  Pd  in  the  diffusion  region  equals  Pt  in 
the  well-stirred  region,  Pd  in  the  diffusion  region  remains  unchanged  and  all  changes 
associated  with  incorporation  of  the  added  volume  can  be  limited  to  the  development  of 
perfusion  heterogeneity  extending  continuously  from  the  original  diffusion  region  into 
only  the  added  volume,  although  more  widespread  changes  in  Ps  and  F  are  possible. 
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Similarly,  adjacent  bubbles  exchange  diffusion  region  volume  if  the  flux  in  the  diffusion 
region  of  one  bubble  abutting  the  diffusion  region  of  an  adjacent  bubble  changes 
direction  from  one  bubble  to  the  other.  Sink  pressure  changes  throughout  the  diffusion 
regions  of  both  bubbles  are  averted  only  if  the  spatial  average  gas  tensions  in  the  two 
diffusion  regions  are  the  same  and  unchanged  by  the  volume  exchange  per  se. 

We  therefore  have  that  (Pd)allk  =Pt  =Pt  throughout  the  lifetimes  of  all  bubbles  in  the 
tissue.  This  relationship  must  be  maintained  as  the  bubbles  evolve  by  properties  of  the 
deviations  functions  Fk(r,  0,  cj)),  which  are  generally  functions  of  time  and  can  be 
different  for  each  bubble. 

Thus  assuming  that  s(t)  =  0  as  in  the  3RUT  model,  Eq.  (A125)  becomes 


1  N  j 

1  IAp.v,).. 


a'tVt  dt 


k=l 


(A126) 


Eq.  (A120),  one  for  each  bubble,  together  with  Eq.  (A126)  determine  the  temporal 
course  of  bubble  growth  and  resolution. 

Note  that  the  accommodation  of  multiple  bubbles  distinguishes  Eq.  (A126)  from  its 
3RUT  model  analog,  Eq.  (A98).  In  the  3RUT  model,  the  constant  diffusion  region 
volume  and  zero  flux  at  the  outer  boundary  of  this  volume  preclude  any  interaction 
between  gas  tensions  in  the  diffusion  region  and  the  well-stirred  region,  obviating  any 
interaction  between  different  bubbles  in  a  given  tissue  volume.  In  the  3RUT-MB  model, 
however,  such  interactions  are  forced  to  occur  by  the  implicit  changes  of  the  diffusion 
region  volumes  around  the  bubbles  as  the  bubbles  evolve. 


Bubble  Number  Density  Limit 

The  number  of  bubbles  that  can  exist  in  a  given  volume  of  tissue  at  any  given  time  is 
limited  by  the  finite  amount  of  gas  in  the  tissue.  An  expression  for  the  bubble  number 
density  limit  is  obtained  by  starting  with  the  3RUT-MB  analog  of  Eq.  (A100)  for  a  given 
bubble,  noting  that  the  flux  f0  at  the  outer  boundary  of  the  bubble  diffusion  region  is 
zero: 


abQt  f  [Ps  - P(z> t)]dv  =  fi0=^-fj.  (A1 27) 

Vd  s 

The  following  3RUT-MB  analog  of  Eq.  (A102)  for  the  gas  content  of  the  diffusion  region 
Ud  is  then  obtained  from  Eq.  (A127)  by  evaluating  the  first  term  of  the  integral  after 
substituting  for  Ps  from  Eq.  (A107),  using  the  relationship  i=at/abQt  ,  and  rearranging 
terms: 
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(A128) 


Ud  =at  JP(z,t)dv  =  atVd(pd  +F)-x^fi, 

vd  s 

where  F=—  f  F(p,0,<|))dv  is  the  spatial  average  of  F(r,  0,  (j>),  and  (pd  +f)  is  the  average 

Vh  ■* 
d  vd 

sink  pressure  Ps.  As  in  the  case  of  Eq.  (A121),  the  validity  of  Eq.  (A128)  does  not 
depend  on  the  shape  of  the  diffusion  region. 

As  for  the  single  bubble  in  the  3RUT  model,  but  in  this  case  for  the  kth  bubble,  the 
average  gas  tension  Pd  in  the  diffusion  region  is  obtained  from  the  expression  for  the 
gas  content  of  the  diffusion  region,  Ud=atPdVd. 

Expressing  fj  as  the  rate  of  change  of 
bubble  gas  content,  we  obtain  the  following  from  Eq.  (A128): 

Fk(Vd)k=4^^(pivi)k-  (A129) 

at  Ds  at 

The  implications  of  Eq.  (A129)  for  the  behavior  of  Fk  during  the  lifetime  of  the  kth  bubble 
are  the  same  as  those  described  following  Eq.  (A104)  for  F  in  the  diffusion  region 
around  the  single  bubble  in  the  3RUT  model:  Fk  remains  positive  during  bubble 
expansion  due  to  increasing  bubble  gas  content,  and  turns  negative  during  bubble 
contraction,  eventually  reaching  zero  when  the  bubble  has  completely  resolved. 

Eq.  (A129)  also  specifies  how  the  diffusion  region  volume  (Vd)k  must  change  with  Fk 
during  bubble  growth  and  resolution.  Determination  of  such  changes  requires 
specification  of  the  function  Fk,  but  such  determination  is  not  necessary  to  solve  the  gas 
and  bubble  dynamics  equations  in  this  model.  It  is  important  to  note,  however,  that  such 
changes  can  occur  within  the  scope  of  the  model,  allowing  adjoining  bubbles  to 
exchange  diffusion  region  space  and  dynamically  compete  for  the  available  gas  and 
perfusion. 

It  follows  from  Eq.  (A127)  that  the  average  sink  pressure  ps  =(Pd  +Fk)  must  remain  non¬ 
zero  during  the  lifetimes  of  all  bubbles:  Ps  must  exceed  the  average  local  gas  tension 
during  bubble  expansion  and  must  be  less  than  the  average  local  gas  tension  during 
bubble  contraction.  We  assume  that  (Vd)k  is  at  least  as  large  as  the  minimum  required  to 
meet  this  condition.  As  multiple  bubbles  compete  for  the  available  gas,  it  is  necessary  to 
ensure  that  the  bubble  number  density  does  not  exceed  the  maximum  that  can  be 
accommodated. 


N  _  _  _  N  _ 

If  (Pd  +Fk)>0  for  each  k,  then  ]T(Pd  +Fk)(Vd)k  =PdVt  +^Fk(Vd)k  >0.  Using  Eq.  (A129), 

k=l  k=l 

this  inequality  becomes 
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(A130) 


Db 

Ds 


The  second  term  of  the  inequality  in  Eq.  (A140)  involves  all  bubbles  present  and  cannot 
exceed  the  total  available  tissue  gas  content  PdVt  at  any  instant,  thus  limiting  the  bubble 
number  density.  Note  that  the  rate  of  change  in  bubble  gas  content  for  each  bubble 
d(PjVj)k/dt  can  be  evaluated  using  Eq.  (A1)  with  the  gradient  gi  given  by  Eq.  (A1 19). 


A-43 


APPENDIX  B.  CORRESPONDENCE  BETWEEN  DIFFUSION  AND  PERFUSION 

EQUATIONS 


In  models  of  gas  bubble  dynamics  in  unstirred  tissue,  a  sink  term  in  the  diffusion 
equation  accounts  for  the  influence  of  tissue-blood  gas  exchange  on  the  diffusion  of  gas 
between  tissue  and  bubble.  The  sink  term  accounts  for  gas  lost  or  gained  by  the 
perfusing  blood,  as  the  gas  flux  term  in  the  perfusion  equation  accounts  for  gas  lost  or 
gained  by  the  bubble.  We  show  here  that  because  both  equations  are  based  on  mass 
balance,  the  perfusion  equation  is  simply  an  integrated  version  of  the  diffusion  equation. 
The  correspondence  between  the  diffusion  and  perfusion  equations  helps  identify  the 
proper  representation  of  the  sink  term  and  ensure  consistency  of  the  model  equations. 

General  Forms  of  the  Equations 

Diffusion  Equation 

The  diffusion  equation  that  governs  the  tissue  gas  tension  in  an  asymmetric  tissue- 
bubble  system  with  a  bubble  centered  at  z=0  is  given  generally  by 

at^^-atDtV2P-abQt[l  +  q(z)][P(z,t)-{Pd(t)  +  F(z,t)}],  (B1) 

dt 

where  z  represents  the  radial,  polar,  and  azimuth  coordinates  (r,  0,  §)  of  an  arbitrary 
point  in  the  diffusion  region  around  the  bubble,  and 

P(z,t)  =  Tissue  gas  pressure  at  point  z  and  time  t, 

at  =  Solubility  of  gas  in  tissue, 

ab=  Solubility  of  gas  in  blood, 

Dt  =  Diffusivity  of  gas  in  tissue, 

Qt=  Blood  flow  rate  into  and  out  of  tissue, 

q(z)  =  Fraction  of  (non-uniform)  blood  flow  at  point  z,  and 

Pd (t)  +  F(z,  t)  =  Ps (z,  t)  =  Sink  pressure  with  component  Pd(t)  that  is  independent  of 

z  and  component  F(z,t)  that  varies  with  z. 

Note  that  the  sink  pressure  defined  here  is  more  general  than  the  sink  pressure  in  the 
3RUT  and  3RUT-MB  models,  where  Pd(t)  is  taken  to  be  the  spatial  average  tissue  gas 

tension  Pd(t)  [see  Eqs.  (A58)  and  (A107)].  The  present  general  definition  reflects  that 
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the  z-independent  pd(t)  may  be  chosen  conveniently  without  any  constraint  to  simplify 
the  model  equations  so  long  as  the  z-dependent  F(z,  t)  is  non-zero. 

The  term  Qt[i  +  q(z)]  accounts  for  variations  of  blood  flow  over  time  and  space.  To 
maintain  the  same  blood  flow  rate  into  and  out  of  the  tissue  with  no  net  accumulation  or 
depletion  of  blood  in  the  tissue,  we  require  the  spatial  variations  represented  by  q(z)  to 
sum  to  zero  over  the  diffusion  region  volume,  Vd;  that  is 

Jq(z)dv  =  0.  (B2) 

vd 


Eq.  (B2)  is  assumed  to  hold  even  in  situations  in  which  blood  flow  varies  with  time,  such 
as  during  exercise,  and  that  such  temporal  changes  in  blood  flow  do  not  alter  the  spatial 
distribution  of  blood  flow.  Dividing  both  sides  of  Eq.  (B1 )  by  at  yields 

=  DtV2P-— [l  +  q(z)]  [P(z,  t)  —  {Pd  (t)  +  F(z,  t)}] ,  (B3) 

at  t 

where  x  =  at/abQt  is  the  blood-tissue  gas  exchange  time  constant. 

Integration  of  Eq.  (B3)  over  the  diffusion  region  volume  yields 

dv  =  J  (d t  V 2p)dv  J  [l  +  q(z ;)] [P(z,  t)  -  {Pd  (t)  +  F(z,  t)}]dv , 


i.e., 


v. 


1  rdP(z,t) 

Vd  l  dt 


—  f (atDt V2p)dv  -if  [P(z, t)  -  Pd (t)]dv 

CX.  f  T 


(B4) 


i  f  q(z)P(z,  t)dv  +  fq(z)dv  +  i  fF(z,t)dv  +  i  f  q(z)F(z,t)dv 

T  J  t-J  t  *  t  ^ 


_  J  j* 

Defining  Pd  = — d  P(t)dv=  Average  tissue  gas  tension  in  the  diffusion  region, 

VH  J 
d  vd 

f0  =  Gas  flux  into  the  tissue  diffusion  region  at  the  outer  boundary  of  the 
diffusion  region, 

fi  =  Gas  flux  into  the  gas  bubble  from  the  tissue  diffusion  region, 
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p  =_L 
q  vd 


j"  q(z)P(z,  t)dv ,  F  =  -j—  |  F(z,  t)dv  ,  Fq  =  -J-  Jq(z)F(z,t)dv  , 

VH  Vd  Vd  Vd  VH 


noting  that  J q(z)dv  =  0  per  Eq.  (B2),  and  evaluating  the  integrals  in  the  first  two 


terms  on  the  right  side,  Equation  (B4)  becomes 


dPd(t)  f0-fj  Vdr- 


[Pd  (0  -  Pd  W]  -  —  Pq  +  —  F  +  —  Fq 


dt 


af 


(B5) 


All  models  considered  in  present  work  postulate  a  zero  gas  flux  at  the  outer  boundary  of 
the  diffusion  region,  i.e. ,  that  f0  =  0.  Thus,  dividing  by  Vd,  using  the  definition  of  the 
blood-tissue  gas  exchange  time  constant,  and  rearranging  terms  yields 


dPd(t)_  Pd(t)-Pd(t) 


dt 


atVd 


Fq  F  Fq 
- -  +  —  +  — 


(B6) 


Perfusion  Equation 

Referring  to  Eq.  (A121),  the  mass  balance  equation  for  blood-tissue  gas  exchange  in 
the  diffusion  region  of  volume  Vd  around  a  bubble  is  given  generally  by 

ai  JP(z,t)dv  =  ab  jQt[(l  +  q(z,t)][Pa  -Pv(z,t)]dv-f1 ,  (B7) 

dt  Vd  vd 

where  we  have  replaced  Qn(z,t)  and  P(z,t)  on  the  right  side  of  Eq.  (A121)  by 
Qt[(l  +  q(z,t)]  and  Pv(z,t),  respectively.  Pv(z,t)  represents  the  venous  gas  tension,  which 
is  allowed  to  vary  here  with  position  in  the  diffusion  region  without  equilibration  with  the 
tissue  gas  tension.  Integrals  on  both  sides  of  Eq.  (B7)  lead  to: 


atVd 


dPd(t) 

dt 


[Pa(t)-Pv(zA)]dv  +  abQt  J  q(z)[Pa(t)-Pv(z,t)]dv-fi 

vd 


=  abQt  J  [PaW-PfoOjdv  +  abQt  J  [P(z,  t)  -  Pv  (z,  t)]dv  + 

Vd  Vd 

(B8) 

«bQt  j  q(z)[Pa(t)ti,pv(z»t)]dv-fi 

Vd 
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Evaluating  the  integral  in  the  first  term  on  the  right  side  of  Eq.  (B8)  and  dividing  both 
sides  by  atvd,  we  get 


J  d[P(z,  t)  -  Pv  (z,  t)]dv  +  J  q(z)  [Pa  (t)  -  Pv  (z,  t)]dv  -  — .  (B9) 

«  T'  'r'/d  yd  tvd  ^  at  vd 


Noting  that  J  q(z)Pa (t)dv  =  Pa (t)  J  q(z)dv  =  0  per  Eq.  (B2),  Equation  (B9)  becomes 

vd  vd 


dPd(t)  _  Pa(t)-Pd(t) _ f,  (  Pd(t)  Pv  Pyg 

dt  T  atVd  TXT 


(BIO) 


where  Pv  =-f-  f  Pv(z,t)dv,  and  p  =-f-  f  q(z)Pv(z,t)dv. 

Vj  V,  V1  i 

Eq.  (BIO)  is  derived  by  considering  that  the  conditions  governing  gas  exchange 
between  blood  and  tissue  are  identical  to  those  presumed  to  integrate  the  diffusion 
equation  and  obtain  Eq.  (B6).  Thus,  Eq.  (BIO)  is  redundant  as  it  is  simply  the  integrated 
version  of  the  diffusion  equation,  which  therefore  accounts  for  all  mass  balance  of  gas 
in  the  tissue.  In  fact,  Eq.  (BIO)  is  identical  to  Eq.  (B6)  if  the  z-independent  sink  pressure 
component  equals  the  arterial  gas  tension,  i.e.,  if  Pd(t)  =  Pa(t),  and  F(z,  t)  and  q(z)  are 
specified  so  that  [f  +  Fq  -Pqj=[pd(t)-  Pv  -  Pvqj. 


Homogeneous  Perfusion 

In  the  case  of  homogeneous  perfusion,  q(z)  =  0  andFa(z,t)  =  o,  and  consequently, 

Pd(t)  =  Ps(z,t)  at  all  z  and  Pq ,  F ,  and  Fq  are  all  zero.  The  presumption  of  homogeneous 

perfusion  implies  that  the  sink  pressure  Ps(t)  must  equal  the  arterial  gas  tension  Pa(t)  at 
every  point  in  the  tissue.  Eq.  (B6)  thus  reduces  to 

dPd(t)  _  Pa  (t)-Pd(t) _ fi__  (B11) 

dt  x  atVd 

Similarly,  we  assume  that  J  [P(z,  t)  -  Pv (z,  t)]dv  =  0 ,  which  implies  that,  on  the  average, 

vd 

the  difference  between  tissue  and  venous  gas  tensions  is  zero.  This  assumption  is  less 
stringent  than  the  assumption  normally  made,  i.e.,  that  the  tissue  is  in  equilibrium  with 

P  ft)  P 

venous  blood  and  P(z,t)  =  Pv(z,t)  at  all  points  in  the  tissue.  Thus  — ^  =  o.  Also, 

X  X 
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because  q(z)  =  0,  Pvq=0  and  Eq.  (BIO)  reduces  to  Eq.  (B1 1)  as  well.  The  requirements 
for  Eqs.  (B6)  and  (BIO)  to  be  identical  are  therefore  satisfied  under  these  conditions. 

With  the  deviation  F(z,  t)  being  zero  for  a  tissue  with  homogeneous  perfusion,  the 
diffusion  equation,  Eq.  (B3),  reduces  to  o 

=  Dtv2P  -  — [p(z,  t)-Ps  (t)],  (B12) 

Ot  T 

which  has  no  general  analytic  solution. 

The  volume  of  tissue  is  absent  in  Eq.  (B12),  as  would  be  expected  with  mass  balance 
being  the  same  in  every  unit  volume  of  tissue.  However,  tissue  volume  appears 
explicitly  in  the  solution  of  Eq.  (B12)  based  on  the  quasi-static  approximation,  under 
which  5P(z,t)/a  on  the  left  side  of  Eq.  (B12)  is  zero.  We  then  have 

DtV2P-  — [P(z,  t)-Ps(t)]  =  0  ,  (B13) 


where  the  sink  pressure,  Ps(t)  =  Pa(t),  need  not  be  equal  to  the  venous  tension  Pv(t). 
Integrating  Eq.  (B13)  over  the  diffusion  region  volume,  we  obtain 


l 

at 


f 


atDtV2Pdv  -  —  f 
t  X 


[P(z,t)-Ps(t)]dv  =  0, 


i.e.,  — — — [Pd(t)-Ps(t)]=0-  (B14) 

at  i 

With  zero  flux  at  the  outer  boundary  of  the  diffusion  region,  fo=0  and  Eq.  (B14)  reduces 
to 


?s(t)  Pd(t) - =  0,  (B15) 

x  atVd 

which  is  inconsistent  with  Eq.  (B1 1 )  if  Ps(t)  =  Pa(t),  implying  that  dPd(t)/dt  =  0.  Therefore 
Ps(t)  cannot  equal  Pa(t)  with  a  finite  tissue  volume  under  the  quasi-static  approximation. 
Moreover,  the  boundary  condition  that  dP(z,t)/dr  =  o  at  the  finite  outer  boundary  of  the 
diffusion  region  relates  pd(t) ,  fj,  and  Vd  through  Eqs.  (A96)  and  (A100),  which  are 

dP  l  r° 

combined  with  e(t)  =  0 ,  fo=0,  and  Ds=Db  to  yield  — -  = - [ [pa (t)  -  Ps (t)]47tp2dp  .  This 

dt  TVd  J 

ri 

relationship  cannot  be  simultaneously  satisfied  with  Eq.  (B11)  unless  Ps(t)  =  Pd(t)  and 
vd  =  00  ■  Thus,  the  only  feasible  solution  to  Eq.  (B13)  is  with  an  infinite  diffusion  region 
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volume  and  Ps(t)  =  Pd(t).  This  satisfies  both  Eqs.  (B1 1 )  and  (B15)  as  well  as  the 
boundary  condition  for  Eq.  (B1 3)  at  r  =  oo ,  which  is  dP(z,  t)/dr  =  0  ,  leading  to 
P(z,t)  =  Ps(t)  =  Pd(t) .  In  what  follows  we  show  that  a  finite  tissue  volume  is 
accommodated  if  we  assume  the  perfusion  to  be  heterogeneous,  which  is  more 
realistic. 


Heterogeneous  Perfusion 

In  the  case  of  heterogeneous  perfusion,  q(z)^0  and  F(z,t)*0,  but  Eq.  (B6)  based  on 
the  diffusion  process  must  still  be  identical  to  Eq.  (BIO)  based  on  blood-tissue  gas 
exchange  for  theoretical  consistency.  As  indicated  earlier,  the  requirements  for  this 
identity  are  Pd(t)  =  Pa(t)  and 


F  +  Fq-Pq=Pd-Pv-Pvq.  (B16) 

Eq.  (B16)  implies  (suppressing  time  dependency  for  notational  simplicity) 

{  (F(z)  +  q(z)F(z)  -  q(z)P(z)}dv  =  J  {P(z)  -  Pv(z)  -  q(z)Pv(z)}dv  , 

Vd  Vd 

i.e.,  J  |F(z)  +  q(z)F(z)}dv  =  J  {P(z)  +  q(z)P(z)-Pv(z)-q(z)Pv(z)}dv, 

Vd  Vd 

i.e.,  J  {l  +  q(z)}F(z)dv  =  J  {l  +  q(z)}{P(z) - Pv (z)}dv  .  (B17) 

Vd  Vd 


Thus,  the  term  F(z)  in  the  diffusion  equation  accounts  for  differences  between  tissue 
and  venous  gas  tensions  caused  by  perfusion  heterogeneity.  If  this  heterogeneity  is 
ignored,  q(z)=0,  F(z)=0,  and  Eq.  (B17)  reduces  to  J  |P(z) - Pv (z)}dv  =  0 ,  which  implies  that 

Vd 

on  the  average,  the  difference  between  tissue  and  venous  gas  tensions  is  zero,  and 
proves  the  assumption  made  earlier  to  reduce  Eq.  (BIO)  to  Eq.  (B1 1 ).  If  only  q(z)  is 

zero,  Eq.  (B17)  still  applies  but  the  term  abQt  J  q(z)[Pa(t)-Pv(z,t)]dv  in  Eq.  (B8) 

vd 

vanishes,  implying  that  the  net  gas  content  carried  by  blood  into  or  out  of  the  diffusion 
region  by  the  heterogeneous  perfusion  components  is  zero. 

In  summary,  we  have 

d¥(z,t)  _ Dtv2p  _l[i  +  q(z)][p  -  {Pa  +F(z,t)}]  for  diffusion  (B18) 

dt  x 
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and 


— -  =  — - - - - —  + - - - 3-  for  perfusion.  (B19) 

dt  t  atVd  t 

Eq.  (B19)  for  perfusion  is  exactly  the  integrated  version  of  Eq.  (B18)  for  diffusion  with 
zero  flux  at  the  outer  boundary  of  the  tissue  diffusion  region  (f0  =  0).  Eq.  (B19),  and 
hence  Eq.  (B18),  remain  consistent  with  mass  balance  between  blood  and  tissue  if  q(z) 
and  F(z)  are  related  to  Pv(z)  as  given  by  Eq.  (B17). 

Analytic  solutions  of  these  equations  are  not  feasible  but  numerical  solutions  can  be 
considerably  streamlined  without  severely  compromising  model  generality  by  simplifying 
the  expressions  for  q(z)  and  F(z,t).  However,  to  ensure  proper  mass  balance  regardless 
of  how  these  expressions  are  formulated,  it  is  imperative  that  the  diffusion  equation 
accounts  for  gas  exchange  between  tissue  and  the  blood  perfusing  the  tissue  and, 
conversely,  that  the  perfusion  equation  accounts  for  gas  exchange  between  tissue  and 
bubbles. 


Sink  Pressure  and  the  Quasi-Static  Approximation 

The  3RUT  and  3RUT-MB  models  are  based  on  the  diffusion  Eq.  (B18)  without  the 
3P(z,t)/St  term  per  the  quasi-static  approximation,  but  the  model  equations  do  not 
explicitly  involve  the  deviations  function  F(z,  t).  Therefore,  the  unspecified  F(z,  t)  may  be 
assumed  to  account  for  the  neglect  of  the  dP(z,t)/<3t  term,  thus  diminishing  the  error  due 
to  the  approximation  in  the  solution. 

Let  F(z, t)  =  Fj (z, t)  +  F2 (z, t) .  Substituting  for  F(z,t),Eq.  (B 18)  becomes 

=  D(V2P  -i[l  +  q(z)][P  -  {Pd  +  F,  (z,  t)  +  F2(z,t)}].  (B20) 

dt  x 

If  we  assign  ap(zT)  =  — [l  +  q(z)]F2 (z, t)  ,  (B21) 

dt  x 

then  Eq.  (B20)  reduces  to 

0  =  DtV2P  -  — [l  +  q(z)][P  -  {Pd  +  F,  (z,  t)}]  -  (B22) 

x 

Eq.  (B22)  replaces  the  unknown  function  F(z,t)  in  Eq.  (B18)  by  another  unknown 
function  Fj(z,t).  More  importantly,  it  eliminates  the  5P(z,t)/5t  term  from  Eq.  (B18)  without 
invoking  the  quasi-static  approximation.  Although  the  partitioning  of  the  deviations 
function  F(z,t)  into  two  components  circumvents  the  quasi-static  approximation,  it 
cannot  completely  eliminate  the  effects  of  the  approximation  on  the  solution.  While  the 
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parameter  X  in  the  3RUT  model  is  a  constant  related  to  the  gas  diffusivity  and  the  blood- 
tissue  gas  exchange  time  constant,  the  parameter  A  in  the  3RUT-MB  model  is  related  to 
the  spatial  average  of  the  deviations  function  and  is  generally  a  function  of  time. 
Therefore,  it  is  desirable  to  treat  A  as  function  of  time  in  the  3RUT-MB  model  in  order  to 
minimize  errors  that  can  arise  from  the  quasi-static  approximation. 
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APPENDIX  C.  EFFECTS  OF  TISSUE  ELASTICITY  ON  BUBBLE  PRESSURE 


The  Young-Laplace  equation  for  the  pressure  in  a  mechanically  stable  gas  bubble  must 
be  modified  to  account  for  the  influence  of  elastic  effects  arising  from  any  rigidity  of  the 
surrounding  medium.  A  term  adopted  from  Gernhardt7  based  on  the  tissue  modulus  of 
compression  is  included  in  Eq.  (A2)  to  account  for  such  influence  on  bubbles  in  tissue. 
Goldman19  correctly  pointed  out  that  this  representation  is  theoretically  inconsistent  and 
proposed  an  alternative  to  Eq.  (A2)  in  which  elastic  effects  on  bubble  pressure  are 
expressed  in  terms  of  the  shear  modulus.  However,  the  analysis  leading  to  the 
alternative  expression  is  not  only  fundamentally  flawed  but  also  neglects  the  steep 
increases  in  shear  modulus  that  can  occur  as  bubble  expansion  causes  approach  to  a 
limiting  stretch  in  hyperelastic  materials  such  as  tissue.  We  here  show  that  Goldman’s 
proposal  is  based  on  an  inappropriate  application  of  infinitesimal  strain  theory  and 
describe  a  correct  alternative  to  Eq.  (A2)  based  on  hyperelastic  theory,  an  alternative 
that  is  roughly  approximated  by  Eq.  (A2)  for  some  systems. 


Infinitesimal  Strain 


We  follow  the  theory  of  infinitesimal  elastic  deformation  presented  by  Landau  and 
Lifshitz20  and  exemplified  without  surface  tension  effects  in  Problem  2  of  their  text.  The 
assumption  that  all  deformations  and  deformation  gradients  are  infinitesimally  small  is  a 
key  feature  of  the  theory  that  allows  linearization  of  the  finite  strain  tensor  by  the  neglect 
of  all  second-order  terms.  With  this  linearization,  the  condition  for  mechanical 
equilibrium  in  a  spherically  symmetric  deformation  considered  in  spherical  coordinates 
with  origin  at  the  center  of  the  material  is 

divu  =  r“2d(r2ur)/ dr  =  constant  =3a,  (Cl) 


with  solution 

ur  =  ar  +  b/r2  ,  (C2) 

where  u  is  the  displacement  vector  of  the  material,  ur  is  the  radial  component  of  u  ,  and 
a  and  b  are  constants  characteristic  of  the  material.  The  principal  components  of  the 
strain  tensor  uik  are  then  given  by 

Uj,.  =  5ur /8r  =  a  -  2b/r3  (C3.a) 


and 


U00  -u<|>c|> 


=  ur  /r  =  a  +  b/r3 . 


(C3.b) 
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Under  the  same  assumption  of  infinitesimal  deformation  and  the  added  condition  that 
the  strains  remain  in  the  linear  range  where  Hook's  law  applies,  the  stress  tensor 
components  are  generally  given  in  terms  of  the  corresponding  strain  tensor  components 
by9 


CT 


ik 


(l  +  x)LU‘k  (1-2  x) 


u,,5 


11  °  ik 


(C4) 


where  E  is  the  Young's  modulus  of  the  material,  x  is  the  Poisson's  ratio  of  the  material, 
and  5ik  is  the  Kronecker  delta. 


We  now  assume  a  spherical  bubble  of  radius  r;  to  be  at  the  center  of  a  spherical  shell  of 
isotropic  elastic  material  with  outer  boundary  at  Ro.  Applying  the  mechanical  equilibrium 
condition  in  Eq.  (Cl ),  the  radial  component  of  the  stress  tensor  associated  with  the 
strain  tensor  is  obtained  by  substituting  Eqs.  (C3.a)  and  (C3.b)  into  Eq.  (C4)  to  yield 

=  3aK -  (4bG / r 3 )  ,  r;<r<R0,  (C5) 

where  K  =  E/3(l-2y)  is  the  modulus  of  compression  of  the  elastomer,  G  =  E/2(l  +  y)  is 
the  small  stress  shear  modulus  of  the  elastomer,  and  the  constants  a  and  b  are 
determined  from  the  boundary  conditions: 

PB  =-CTrr(ri)+2Y/E  (C6) 


and 


P0=-arr(R0),  (C7) 

where  PB  is  the  gas  pressure  in  the  bubble,  an(ri)  is  the  radial  stress  in  the  elastic 
material  at  the  bubble  surface,  P0  is  the  hydrostatic  pressure  acting  on  the  elastic 
material  at  distance  R0,  arr(R0)  is  the  radial  stress  in  the  elastic  material  at 
distance  R0,  and  y  is  the  surface  tension  of  the  bubble-elastomer  interface. 

Goldman19  takes  up  the  problem  at  this  point  in  the  analysis  and,  in  a  departure  from 
other  treatments,  applies  the  mechanical  equilibrium  condition  in  Eq.  (Cl)  with  solution 

given  by  Eq.  (C2)  to  the  gas  in  the  bubble  (u1g)  =  a(g)r +  b(g) /r2 ).  He  then  uses  the 

continuity  of  the  displacement  vector  at  the  bubble  surface, 


9  We  adopt  the  usual  convention  for  representing  tensors:  Whenever  a  suffix  is 
repeated  in  a  given  term  it  is  to  be  given  all  possible  values  r,  0,  and  4>  and  the  terms  are 
to  be  added  for  all. 
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U  rg)  (ri  )  =  U  r  (r;  ) , 


(C8) 


as  another  boundary  condition  to  obtain 

b  =  r;3(a(g)  -  a),  (C9) 

where  u(rg)(rj)  is  the  displacement  vector  in  the  gas  and  ur(r;)  is  the  displacement 
vector  in  the  elastomer,  each  at  the  bubble  surface,  and  a(g)  and  b(g)  (b(g)  =  0  )  are 
constants  characteristic  of  the  gas  in  the  bubble.  Goldman  argues  that  3a(g)  is  the  trace 
of  the  strain  tensor  u[^  of  the  gas  in  the  bubble  given  in  terms  of  the  relative  volumetric 
strain  of  the  gas: 

uf  =(u£>  +u<i)  +  u<*))=Av/v0  =(v-v0)/v0  (CIO. a) 


=  3a(g) . 


(ClO.b) 


where  Eq.  (ClO.b)  follows  from  Eq.  (Cl ).  Goldman  identifies  v0  as  the  volume  of  gas  in 
its  undeformed  state  at  zero  pressure  ( p0  =  o )  and  V  as  the  final  volume  of  the  gas  in 
the  bubble  at  pressure  PB  in  the  deformed  elastomer.  Thus  as  p0  ->  o ,  V0  ->oo  and  from 
Eq.  (CIO) 

a(g)  — >  —1  / 3  ,  (Cl  1 ) 


regardless  of  V  .  This  result,  a(g)  =  -1/3 ,  is  then  used  with  Eqs.  (C5),  (C7),  and  (C9)  to 
obtain  the  following  solution  of  Eq.  (C6): 


PB  = 


1  +  (4G/3K) 
1  +  (4Gv/3K) 


|P0  + 


4Gv^ 


4G 

3 


+  2y  /  r; , 


(C12) 


where  v  =  r;3/Ro  and  P0  is  as  defined  for  Eq.  (C7).  For  an  incompressible  elastomer, 

K  =  oo  and  Eq.  (Cl  2)  reduces  to 

PB=Po+^-(v-l)  +  2y/r].  (C13) 

Because  (v-l)<0  always,  Eq.  (Cl  3)  implies  that  shear  forces  always  exert  negative 
effects  on  the  bubble  pressure  in  an  incompressible  elastomer. 

The  analysis  leading  to  Eqs.  (Cl  2)  and  (Cl  3)  is  flawed  because  Eqs.  (Cl )  and  (Cl  O.a), 
and  hence  Eq.  (ClO.b),  are  applicable  only  to  very  small  deformations,  not  to  the 
deformation  of  a  mass  of  gas  from  an  effectively  infinite  volume  to  the  finite  volume  of  a 
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real  bubble.  The  importance  of  this  constraint  is  illustrated  in  the  derivation  of  Eq. 

(CIO. a)  given  by  Landau  and  Lifshitz.20  For  a  deformation  in  which  dx;  is  the  radius 
vector  joining  two  points  in  a  material  along  the  ith  principal  axis  before  deformation  and 
dxj  is  the  radius  vector  joining  same  two  points  along  the  same  axis  after  deformation, 
dxj  is  given  in  terms  of  the  ith  principal  component  u(i)  of  the  strain  tensor  by 


dx;  = 


1  + 


+  2u(i))  -1 


dx- 


(C14) 


If 


u 


(i) 


«1 


,  all  multiples  of  u(i)  can  be  neglected,  fA/(l  +  2u(1)) -ll  =  u(1) ,  and  Eq.  (C14) 


becomes 


dxj  =  l  +  u(i)  dxj  . 


(C15) 


The  volume  of  the  material  after  the  deformation  is  then  given  by 


dV  =  dV(l  +  Ujj) , 


(01 6) 


from  which  Eq.  (CIO. a)  follows.  However,  all  terms  in  the  strain  tensor,  including 
second-order  terms,  and  absolute  values  of  the  principal  components  of  the  strain 
tensor  approach  their  maximum  possible  values  -  not  negligibly  small  values  -  in  a 
deformation  of  a  mass  of  gas  from  an  effectively  infinite  volume  to  the  finite  volume  of  a 
real  bubble.  Therefore,  neither  Eq.  (Cl)  nor  Eq.  (CIO)  applies  and  no  ready  relationship 
exists  between  the  volumetric  strain  and  the  strain  tensor  for  such  a  deformation. 


Hyperelastic  Strain 


Eq.  (C6)  is  conventionally  solved  using  hyperelastic  theory  in  which  Eq.  (C5)  is  replaced 
by  a  stress-strain  relationship  obtained  from  the  free  energy  of  the  bubble-induced 
deformation  of  the  elastomer,21,22  a  relationship  that  is  independent  of  properties  of  the 
gas  in  the  bubble.  For  a  spherical  bubble  of  radius  r  in  a  homogeneous  isotropic 
incompressible  elastomer  and  boundary  conditions  given  by  Eq.  (C6)  and  Eq.  (C7)  with 
R0  =oo  ,  Zhu,  et  al.,22  present  the  following  solution  for  PB  based  on  the  Gent  free  energy 
function:23 


PB=Po 


+  21  +  2Gp/,'_M1r%_i 

ri  Jl  1-(2A.2  +r4  -  3 )/  J  Iim 


(C17) 
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where  X,  =  ri/ri°  is  the  hoop  or  tangential  stretch  in  the  elastomer  at  the  bubble  surface, 
r°  is  the  radius  of  the  unstretched  bubble,  and  Jnm  is  a  constant  related  to  the  limiting 
stretch  A,lim  by 


'lim 


=  24m+^l-3. 


lim 


(C18) 


The  limiting  stretch  is  the  point  at  which  a  cross-linked  polymeric  elastomer  can  no 
longer  accommodate  increasing  stretch  by  alignment  of  its  constituent  polymer  chains  in 
the  direction  of  the  stress.  The  stress  increases  steeply  as  the  limiting  stretch  is 
approached  and  the  elastomer  fails  as  the  limiting  stretch  is  exceeded. 


With  finite  Jlim,  the  integral  in  Eq.  (Cl  7)  is  evaluated  numerically  to  determine  the 
pressure  in  the  bubble.  For  a  neo-Hookean  material,  jlim  ->oo  ,  and  the  integral  in  Eq. 
(Cl  7)  is  evaluated  to  yield  24 


p  -p  +2UG 

—  *r\  - '  I 


( 

4 

r; 

r; 
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1 

— 

1 

r. 

v  1  / 

v  1  ) 

(C19) 


Eq.  (Cl 7)  or  Eq.  (Cl  9)  can  be  applied  to  each  of  many  bubbles  in  a  given  volume  of 
elastomer  if  each  bubble  can  be  considered  to  be  at  the  center  of  its  own  independent 
spherical  mechanical  domain.  In  such  cases,  bubble  number  densities  are  assumed  to 
be  low  enough  so  that  the  bubbles  do  not  mechanically  interact. 

If  the  bubble  is  unstretched,  r,  =  r° ,  A,  =  l,and  pB  =  p0  +27/^° .  For  bubbles  of  radius 
r  >  r° ,  the  Laplace  pressure  due  to  surface  tension  decreases  monotonically  as  the 

radius  of  the  bubble  increases,  while  the  pressure  due  to  the  elasticity  of  the  elastomer 
is  always  positive  and  monotonically  increasing.  This  latter  feature  of  elastic  effects  is  in 
contrast  to  the  behavior  prescribed  by  the  second  term  on  the  right  in  Eq.  (Cl  3).  As 
illustrated  in  Figure  C.1,  elastic  forces  when  ^/Gr°  is  relatively  small  override  the 
effects  of  surface  tension  effects  at  values  of  X  near  and  smaller  than  1  and  act  to  resist 
dissolution  of  the  bubble  by  reducing  the  bubble  pressure.  As  the  bubble  radius 
increases,  the  elastic  pressure  asymptotically  approaches  5G/2  in  the  neo-Hookean 
material.  In  the  material  with  a  limiting  stretch,  the  elastic  pressure  increases  steeply  as 
the  limiting  stretch  -  and  material  failure  -  is  approached. 
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Figure  C.l.  Comparison  of  bubble  pressures  obtained  with  Eqs.  (Cl 7),  (Cl  9),  and 
(A2):  p0  =  latm,  G  =  1  atm,  r°=  10  pm,  y  =  70  dyne/cm,  jlim=  80  (A,lim=  6.44),  and 
M  =  IxlO7  atm. 


Figure  C.2  illustrates  the  behavior  of  the  bubble  pressure  in  an  elastic  material  with  a 
larger  A./Gr°  ratio  and  finite  limiting  stretch.  In  such  cases,  surface  tension  effects 
predominate  at  low  and  fractional  stretch.  As  the  fractional  stretch  is  further  decreased, 
bubble  pressure  reaches  a  maximum  and  falls  as  elastic  forces  act  to  resist  dissolution 
of  the  bubble. 
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Figure  C.2.  Comparison  of  bubble  pressures  obtained  with  Eqs.  (Cl 7)  and  (A2):  P0  = 
latm,  G  =  0.001  atm,  rj°=  1  pm,  y=  70  dyne/cm,  Jlim=  80  (X,lim=  6.44),  and 
M  =  7x1 09  atm. 


Elastic  effects  in  either  Eq.  (Cl 7)  or  Eq.  (Cl  9)  are  expressed  in  terms  of  the  small 
stress  shear  modulus  G,  not  the  elastic  modulus  M  as  in  Eq.  (A2).  Figure  C.1  illustrates 
that  Eq.  (A2)  provides  a  poor  representation  of  elastic  effects  on  bubble  pressure  in  a 
material  with  small  A,/Gr° .  However,  Figure  C.2  illustrates  that  Eq.  (A2)  can  provide  a 
rough  single-parameter  approximation  of  the  behavior  prescribed  by  Eq.  (Cl 7)  for  a 
material  with  large  A,/Gr°  if  M  is  treated  as  an  empirical  constant. 
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